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1.  PROBLEMS  STUDIED 


This  research  was  motivated  by  the  fact  that  real-time  simulation  of  large  and  complex 
multibody  systems  -such  as  Army  ground  vehicles-  was  not  possible  with  existing  technology 
in  standard  computer  architectures.  This  was  a  challenge  that  needed  to  be  addressed  before 
further  advances  in  mechanical  simulation  with  hardware-in-the-loop  and  man-in-the-loop  (crew 
station),  intelligent  vehicle  control  systems  (IVCS)  and  virtual  prototyping  were  made  possible. 

This  proposed  research  project  was  aimed  at  developing  methods  and  concepts  that  would 
lead  to  faster- than- real-time  simulation  using  standard  computer  platforms  (high-end 
workstations).  One  of  the  concepts  proposed  is  that  of  intelligent  simulation,  which  integrates, 
depending  upon  the  characteristics  of  the  system,  all  the  factors  involved  in  the  simulation 
process.  It  was  envisioned  that  this  concept  would  lead  to  faster  and  more  robust  real-time 
simulators  that  will  strongly  influence  the  missions  of  the  US  Army  TARJDEC  such  as: 
intelligent  vehicle  control  systems  (IVCS),  virtual  prototyping  and  real-time  simulation  stations 
with  soldier-in-the-loop.  This  research  is  in  essence  of  multidisciplinary  nature  encompassing 
engineering  mechanics,  computauonal  science  and  numerical  methods. 

More  specifically  the  following  problems  have  been  addressed: 

•  Develop  robust  and  efficient  numerical  algorithms  in  descriptor  form,  and  hybrid  integration 
methods  for  the  real-time  simulation  of  constrained  mechanical  systems.  These  algorithms 
will  couple  the  dynamic  formulation,  kinematic  constraint  conditions,  with  hybrid 
integration  procedures  to  achieve  methods  that  will  yield  improved  accuracy,  efficiency  and 
robustness. 

•  Establish  the  criteria  and  mechanisms  that  will  allow,  based  on  the  characteristics  of  the 
system  to  be  simulated,  for  the  systematic  choice  and  combination  of  the  modeling, 
formulation  of  the  equations  of  motion,  index  reduction  methods,  numerical  integration  and 
solution  method. 

•  Study  parallelization  schemes  that  will  integrate  the  different  components  that  form  the 
intelligent  simulator  in  order  to  obtain  efficient  solutions  faster  than  real-time. 


2.  SCIENTIFIC  A  CCOMPLISHMENTS  AND  RES  VETS 
2. 1  Formulations  and  benchmark  problems. 

A  library  of  benchmark  problems  has  been  created  that  includes  the  following  systems:  open 
and  close  chain,  medium  and  large  scale  stiff  problems,  systems  with  redundant  constraints, 
high  frequency  response,  and  singular  topologies  (Publications  1  and  3  below).  This  library 
has  served  as  a  testing  ground  for  every  method  previously  formulated  or  developed  in  die 
course  of  this  research  effort. 

Four  state  of  the  art  formulations  have  been  implemented  and  tested  for  simulation:  two 
augmented  Lagrangian  formulations  in  descriptor  form  in  index- 1  and  index-3,  a  new  state- 
space  formulation  suitable  for  stiff  systems,  and  a  new  fuUy-recursive  (order  N)  formulation 
(Publication  3). 

Considering  a  multibody  system  whose  configuration  is  characterized  by  n  coordinates  q 
that  are  interrelated  through  the  m  holonomic  kinematic  constraint  conditions  <J>(q,  0  =  0,  the 
equations  of  motion  become: 

M  q  +  $qT  X  =  Q(q,q) 

where  M  is  the  mass  matrix,  Q  is  the  vector  that  contains  the  external  and  non  conservative 
forces  as  well  as  the  velocity  dependent  inertia  forces,  <I>q  is  the  Jacobian  of  the  constraint 
equations,  and  X  is  the  vector  that  contains  the  Lagrange’s  multipliers.  The  last  two  equations 
constitute  a  set  of  n+m  mixed  differential  algebraic  equations  (DAE)  of  index-3,  with  q  and  X 
as  unknowns.  This  method  is  referred  to  as  the  classical  Lagrangian  formulation. 

Index-1  Augmented  Lagrangian  Formulation  with  Projections 

This  formulation  leads  to  the  following  equations  of  motion: 

Mq  +  +  <l>q  X  =  Q(q,q) 

where  X*  are  the  Lagrange  multipliers.  Introducing  the  expression  ^>  =  C>q  q+  Oq  q  +  O,  in 
equation  the  following  equation  is  obtained 

(M  +  OqOcOq)  q  +  Oq  X  =  Q(q,q)  -  $q  a  (<I>q  q  +  <!>,) 

It  is  important  to  note  that  there  is  a  substantial  difference  between  this  method  and  the 
1  classical  Lagrange's  multiplier  approach  presented  above.  The  leading  matrix  of  the  classical 

1  Lagrangian  method  becomes  singular  in  singular  configurations,  changing  topologies  and  in  the 

(  presence  of  redundant  constraints.  However,  although  the  mass  matrix  M  is  in  general  positive 

>  semi-definite,  the  leading  matrix  of  equation  (M  +  C>Ja<l>q),  is  always  positive  definite,  which 

i  means  that  it  can  always  be  factored,  even  in  singular  positions,  changing  topologies  and/or 
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with  redundant  constraints.  The  details  of  how  the  projections  are  performed  are  given  in  the 
paper  attached  in  the  Appendix. 

Index-3  Augmented  Lagrangian  Formulation  with  Projections 

Publication  3  (see  Appendix)  also  presents  the  index-3  augmented  Lagrangian  method  with 
mass-orthogonal  projections.  Such  formulation  leads  to  the  following  equations  of  motion: 

Mq  +  OqOt  <!>  +  #q  X*  =  Q(q,q) 

where  X*  are  the  Lagrange  multipliers.  The  numerical  implementation  of  this  formulation  also 
leads  to  a  similar  leading  matrix  to  that  of  the  previous  method. 

Modified  State-Space  Formulation 

A  state-space  representation  of  the  equations  of  motion  can  be  obtained  using  the  basis  of  the 
null-space  of  the  Jacobian  of  the  constraints.  Such  basis  is  defined  by  the  matrix  R  which  will 
satisfy  the  following  condition  <J>q  R  =  0.  Accordingly,  the  dependent  velocities  and 
accelerations  can  be  expressed  in  terms  of  the  independent  ones  as  follows: 

q  =  Rz  and  q  =  R  z  +  R  z 

The  end  result  (see  Appendix)  is  the  following  equation: 

RtM  R  i  =  RTQ(q,q)  -RTMRz 

which  constitutes  the  equations  of  motion  in  independent  coordinates.  This  equation  provides 
the  function  evaluation  for  the  numerical  integration  scheme.  This  method  requires  the  solution 
of  the  position  and  velocity  problems  at  each  time  step  to  obtain  the  matrix  R,  q  and  q ,  so  that 
the  accelerations  may  be  obtained.  The  solution  of  the  position  and  velocity  problems,  that  are 
costly  computationally,  have  the  same  meaning  as  that  of  projecting  the  solution  to  obtain  the 
dependent  positions  and  velocities. 

Implicit  integrators  make  successive  calls  to  the  above  equation  as  a  fixed  point  iteration. 
Such  iteration  is  either  generally  inefficient  or  precludes  convergence  for  stiff  systems.  Due  to 
this  reason  a  method  has  been  developed  that  allows  this  method  to  be  used  in  stiff  problems, 
and  is  described  in  detail  in  the  publication  3  (see  Appendix). 

Fully-Recursive  Formulation 

A  topological,  fully- recursive  algorithm  of  order  (9(N)  has  been  developed.  This  method  is 
based  and  improves  on  the  articulated  inertia  method  by  which  accelerations  are  obtained  in  a 
recursive  manner.  The  details  are  presented  in  publication  6,  reproduced  in  the  Appendix. 


Comparative  Results 

Comparative  results,  illustrated  in  Tables  I  and  II,  lead  to  the  following  conclusions  (this  is  the 

first  time  a  comparison  has  been  made  among  the  most  suitable  methods  for  multibody 

dynamics  simulation): 

•  The  index-3  formulationis  very  efficient  but  it  leads  to  serious  ill-conditioning  at  small  time 
steps,  starting  at  about  10  5,  when  the  penalty  parameter  becomes  problem  dependent. 

•  The  index- 1  formulation  with  projections  is  the  most  robust  of  all  the  methods,  but  it  may 
not  converge  for  large  time  steps  where  the  index-3  does.  Contrary  to  index-3  it  becomes 
more  accurate  as  the  time  step  decreases  and  it  is  not  affected  by  ill-conditioning.  It  is  also 
insensitive  to  the  penalty  parameter  which  remained  constant  in  all  the  simulations. 

•  The  new  state-space  method  works  well  with  stiffness,  shows  robust  behavior  and  good 
energy  preservation,  as  well  as  good  convergence  at  all  time  steps.  However,  it  runs 
considerably  slower  than  the  index- 1  and  index-3  methods  in  all  the  cases  tested. 

•  The  fully-recursive  (articulated  inertia)  formulation  behaves  in  a  similar  way  to  that  of  the 
previous  one,  except  in  what  refers  to  speed  and  stiffness.  This  method  is  much  faster  than 
the  state-space.  Also,  its  improves  as  the  size  of  the  problem  increases.  For  large  non-stiff 
problems  the  method  becomes  the  most  competitive. 

•  The  following  tables  summarize  the  findings  of  this  paper.  Each  method  is  given  a  grade 
depending  on  its  performance  under  different  situations.  The  letter  "A”  means  outstanding, 
"B"  good,  "C"  fair,  "D"  poor  and  "F"  signifies  a  failure. 
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Table  II.  Performance  for  systems  under  different  conditions. 


•  Summarizing: 

-  A  multi-index- 1 -index-3  formulation  with  projections  is  the  best  candidate  for  an  efficient 
and  robust  general  purpose  real-time  simulator. 

-  The  fully-recursive  formulation  may  prove  suitable  for  very  large  problems  as  a  special 
purpose  real-time  tool. 


2.2  Modeling  and  performance  of  the  Fully  Recursive,  Augmented 
Lagrangian  and  Classical  Lagrangian  methods. 

The  concept  of  intelligent  simulation  involves  the  coupling  of  the  formulations  with  the 
modeling  and  solution  aspects.  A  key  conclusion  of  this  research  is  the  importance  that  the 
modeling  aspects  has  in  the  overall  efficiency  of  the  simulation  process.  Consequently,  a 
description  is  made  on  some  of  the  modeling  aspects,  as  well  as  a  report  of  the  final  results 
obtained  with  each  of  the  methods  including  numerical  stiffness  in  sequential  as  well  as  in  a 
parallel  processing  environments. 

The  example  used  in  all  the  recent  work  is  the  military  4x4  Bombardier  litis  vehicle  whose 
model  has  been  reported  before.  As  previously  described  the  car  features  four  identical 
suspensions  (see  Figure  1)  composed  of  a  shock  absorber,  which  provides  a  damping  force  Fq 
and  an  elastic  force  Fg  (see  Appendix).  A  leaf  spring,  modeled  as  a  linear  spring  of  35,900 
N/m.  Moreover,  a  rubber  bumper  contacts  the  wheel  at  the  leaf  spring  connection  point  after  a 
vertical  motion,  from  the  equilibrium  position  of  70  mm.  The  bumper  stiffness  is  assumed  to 
be  a  constant  value  of  107  N/m,  which  is  engaged  only  after  the  clearance  is  used  up.  This  last 
element  is  considered  in  some  simulations,  thus  adding  a  much  higher  stiffness  to  the  problem. 
Its  presence  will  be  indicated.  A  tire  is  also  included  whose  radial  stiffness  is  460,000  N/m. 
The  side  force  and  aligning  torque  are  considered  by  application  of  the  tire  Calspan  modci. 
Finally,  there  is  a  steering  system,  which  is  kinematically  guided  during  the  maneuvers.  The 
total  mass  of  the  vehicle  is  about  1500  Kg. 


MODELING  FOR  LAGRAN GIAN  FORMULATIONS 

For  analysis  with  Lagrangian  formulations,  the  litis  vehicle  has  been  modeled  using  fully 
Cartesian  dependent  coordinates,  also  called  natural  coordinates  (see  Appendix).  For  three- 
dimensional  multibody  systems,  these  coordinates  describe  the  position  of  the  whole  multibody 
system  by  means  of  the  Cartesian  coordinates  of  basic  points  and  the  Cartesian  components  of 
unit  vectors,  all  distributed  throughout  the  elements.  Each  body  of  the  system  should  have  a 
sufficient  number  of  points  and  vectors  linked  to  it,  so  that  its  motion  is  completely  defined. 
Figure  1  illustrates  how  the  chassis,  suspension  and  steering  system  of  the  vehicle  have  been 
modeled. 


Figure  I.  Modeling  of  the  litis  vehicle  in  fully  Cartesian  coordinates. 


As  may  be  seen  in  Figure  1,  basic  points,  unit  vectors  and  relative  coordinates  of  distance,  have 
been  used,  the  total  number  of  variables  being  168.  The  system  has  1 1  degrees  of  freedom:  six 
free-body  motion  degrees  of  freedom  for  the  chassis,  one  for  each  suspension  and  one  more  for 
the  steering  system.  In  consequence,  the  variables  of  the  problem  are  related  through  168- 
11=157  constraint  equations.  In  practice,  only  10  from  the  11  degrees  of  freedom  are  true 
dynamic  degrees  of  freedom  because,  as  said  before,  the  steering  motion  is  kinematically 
guided.  This  guidance  is  introduced  through  a  new  constraint  equadon,  so  that  the  final  number 
of  constraints  is  158. 

The  modeling  problem  using  this  type  of  coordinates  becomes  straight- forward  and  it  may 
be  used  in  general  purpose  simulators,  in  fact  its  use  is  very  similar  to  that  of  the  finite  element 
method  in  structural  mechanics. 


MODELING  FOR  RECURSIVE  FORMULATION 

Recursive  formulations  require  the  use  of  relative  coordinates.  They  define  the  position  of  each 
body  in  relation  to  the  previous  one  in  the  kinematic  chain,  by  using  the  parameters  or 
coordinates  corresponding  to  the  relative  degrees  of  freedom  allowed  by  the  joint  linking  those 
elements.  Figure  2  illustrates  the  way  the  chassis,  suspension  and  steering  system  of  the  litis 
vehicle  have  been  modeled  usmg  this  kind  of  coordinates. 


Figure  2.  Modeling  of  the  litis  vehicle  in  relative  coordinates. 

In  this  case,  the  number  of  variables  is  43,  distributed  as  follows:  6  for  the  chassis  (three 
translations  and  three  rotations),  9  for  each  suspension  (all  of  them  rotations  at  the  different 
pairs  as  shown,  in  Figure  2).  and  I  distance  for  the  steering  system.  As  this  last  variable  is 
kinematically  guided,  the  total  number  of  variables  that  are  integrated  is  42.  Figure  2  shows 
where  the  closed-loops  have  been  cut.  Two  cuts  are  performed  for  each  suspension.  One  of 
them  needs  3+3=6  constraints  to  establish  identity  between  points  and  directions  at  the 
rotational  joint  removed  (actually,  only  two  constraints  are  needed  to  specify  the  direction  but 
three  makes  it  easier),  while  the  other  only  introduces  3  more  constraints  as  this  time  the  joint 
removed  is  spherical.  Therefore,  the  total  number  of  constraints  is  36,  with  just  32  independent. 

The  modeling  problem  using  this  type  of  coordinates  is  quite  involved  and  requires  a 
substantial  effort  particularly  in  the  existence  of  close  loops.  Its  use  is  not  recommended  for 
general  purpose  simulators,  but  for  special  ones  that  would  allow  for  a  time  consuming  effort  of 
implementation. 


DESCRIPTION  OF  THE  SIMULATIONS 


The  first  simulation,  called  SIMUL1,  consists  of  10  seconds  of  motion  with  the  vehicle  going 
up  an  inclined  ramp  and  the  down  a  series  of  stairs  at  a  speed  of  5  m/s,  as  may  be  seen  in 
Figure  3.  The  simulation  leads  to  the  strong  motion  with  accelerations  of  up  to  5g. 


Figure  3.  The  litis  vehicle  performing  SIMUL1. 


During  the  second  simulation,  SIMUL2,  that  also  lasts  for  10  seconds  at  a  speed  of  5  m/s,  the 
vehicle  goes  off  a  ramp  and  hits  the  ground  as  illustrated  in  Figure  4. 


Figure  4.  The  litis  vehicle  performing  SIMUL2 


In  the  third  simulation,  SIMUL3,  the  litis  vehicle  is  driven  through  a  series  of  columns  in  the 
slalom  maneuver  shown  in  Figure  5.  The  simulation  lasts  for  10  seconds  at  a  vehicle  speed  of  5 

m/s. 


Figure  5.  The  litis  vehicle  performing  SIMUL3. 


RESULTS  ON  ONE-PROCESSOR  MACHINES.  SEQUENTIAL  PROCESSING. 

The  three  simulations  described  above  were  carried  out  on  a  SGI  Indigo2  IMPACT  with  one 
processor  R4400SC  @  200  MHz  and  2  Mb  of  secondary  cache  memory,  using  the  three 
formulations:  Augmented  Lagrangian,  Classical  Lagrangian  and  Fully  Recursive.  Table  I 
features  two  columns  for  each  formulation:  the  first  one  gives  the  smallest  CPU  time  achieved 
with  the  corresponding  method,  while  the  second  one  shows  the  time  step  size  which  led  to  the 
best  CPU  time. 

Table  III  illustrates  the  fact  that  the  augmented  Lagrangian  formulation  (with  a  penalty 
factor  of  109)  is  clearly  superior  to  the  classical  one  because  it  performs  notably  faster  (between 
three  and  four  times)  and  achieves  convergence  for  greater  time  step  sizes.  On  the  other  hand, 
the  fully-recursive  formulation  shows  to  be  more  efficient  than  the  augmented  Lagrangian 
method  although  it  needs  to  take  a  smaller  time  step  size  to  converge. 

Therefore,  when  the  simulation  difficulty  is  not  very  high  and  the  augmented  Lagrangian 
method  can  work  with  large  time  step  sizes  (SIMUL2),  it  performs  at  the  same  level  that  the 
fully-recursive  formulation.  On  the  contrary,  when  the  maneuver  becomes  sharper  (SIMUL1) 
the  augmented  Lagrangian  method  must  keep  on  small  step  sizes  and  the  fully-recursive 
formulation  takes  a  fair  advantage. 


Augmented  Lagrangian 

Classical  Lagrangian 

Fully-Recursive 

CPU  (s) 

At  (s) 

CPU  (s) 

At  (s) 

CPU  (s) 

At  (s) 

SIMUL1 

34 

0.0175 

139 

0.01 

20 

0.0075 

SIMUL2 

18 

0.03 

58 

0.025 

14 

0.0075 

SIMUL3 

27 

0.025 

99 

0.01 

9 

0.0075 

Table  III.  Efficiency  on  a  single-processor  machine. 


RESULTS  ON  PARALLEL  MACHINES 

To  find  out  the  effect  of  the  parallelization  over  the  three  formulations  under  study,  simulation 
SIMUL1  has  been  executed  on  a  SUN  Sparc  Station  20  with  4  processors  HyperSparc 
(RT625)  @  100  MHz  with  256  Kb  of  cache  memory.  Table  IV  shows  the  reduction  in  CPU 
time  experienced  by  each  method  as  the  number  of  processors  increases. 


Augmented 

Classical 

Fully- 

Lagrangian 

Lagrangian 

Recursive 

1 

0.0 

0.0 

0.0 

2 

22.7 

6.5 

0.0 

3 

29.6 

9.0 

0.0 

4 

29.6 

9.0 

0.0 

Table  IV.  Reduction  of  CPU  time  in  %  due  to  parallelization. 


It  may  be  seen  from  Table  IV  that  the  method  that  takes  most  advantage  of  the 
parallelization,  up  to  30%  savings,  is  the  augmented  Lagrangian,  since  it  uses  it  during  the 
matrix  formation  and  solution.  The  classical  Lagrangian  only  takes  advantage  of  the 
parallelization  at  the  time  of  solving  the  equations,  achieving  up  to  9%  reduction.  On  the  other 
hand,  the  fully-recursive  formulation  does  not  profit  at  all  from  the  parallel  computing,  thus 
indicating  to  be  a  method  suitable  for  sequential  processing.  It  must  be  noticed  that,  in  all  cases, 
the  fourth  processor  does  not  contribute  to  make  any  difference  in  CPU  time. 

The  simulations  SIMUL1  and  SIMUL2  are  also  performed  in  an  ALPHA  AXP  4000  with 
2  processors  Alpha  5  @  466  MHz  with  4Mb  of  cache  memory.  The  performance  of  the 
augmented  Lagrangian  method  is  shown  in  Table  V.  Real-time  performance  is  achieved  for 
SIMUL2,  thus  corroborating  the  fact  that  real-time  simulation  can  be  achieved  in  current 
standard  workstations  using  the  augmented  Lagrangian  formulation. 


Number  of  Processors 

SIMUL1 

SIMUL2 

2 

12  seconds 

9  seconds 

Table  V.  CPU  times  taken  by  the  augmented  Lagrangian  method  in  the  ALPHA  AXP. 


INFLUENCE  OF  NUMERICAL  STIFFNESS 

Simulation  SIMUL1  is  carried  out  again  using  the  three  formulations.  However,  this  time  the 
rubber  bumpers  of  stiffness  107N/m  mentioned  above,  are  activated  in  the  four  suspensions  of 
the  model,  so  that  the  behavior  of  the  algorithms  under  conditions  of  high  numerical  stiffness 
may  be  observed.  This  simulation  is  now  called  SIMUL1S.  Table  VI  compares  the  results 
obtained  by  the  three  methods  when  performing  SIMUL1  and  SIMUL1S. 


.  _  - 

Augmented  Lagrangian 

Classical  Lagrangian 

Fully-Recursive 

CPU  (s) 

At  (s) 

CPU  (s) 

At  (s) 

CPU  (s) 

At  (s) 

SEMUL1 

34 

0.0175 

139 

0.01 

20 

0.0075 

SIMUL1S 

78 

0.005 

435 

0.0025 

34 

0.0025 

Table  VI.  Influence  of  numerical  stiffness. 


Table  VI  shows  that  numerical  stiffness  makes  the  Lagrangian  methods  take  a  small  time 
step  size,  thus  loosing  performance  with  respect  to  the  fully-recursive  formulation.  On  the  other 
hand,  the  efficiency  ratio  between  the  augmented  and  classical  Lagrangian  methods  does  not 
seem  to  change.  However,  the  penalty  factor  of  the  augmented  Lagrangian  procedure  needs  to 
be  increased  to  1010to  achieve  convergence  and,  at  lower  time  steps,  as  for  instance  1(T3 
seconds,  this  formulation  starts  suffering  from  its  characteristic  ill-conditioning  effect.  When 
numerical  stiffness  remains  low,  this  problem  arises  at  a  step  size  of  10-4  or  10-5  seconds, 
but  the  presence  of  high  numerical  stiffness  makes  it  to  appear  sooner. 

CONCLUSIONS 

This  is  the  first  time  that  a  thorough  comparison  of  different  leading  methods  of  dynamic 
simulation  has  been  carried  out  for  large  multibody  systems.  The  conclusions  that  have  been 
reached  may  be  summarized  as  follows: 

«•  The  modeling  process  required  by  the  fully-recursive  method  leads  to  a  relative  coordinate 
formulation  that  is  suitable  for  open  chain  configurations,  but  presents  serious  complexities 
and  difficulties  for  the  case  of  closed  kinematic  chains.  The  modeling  for  the  Lagrangian 
formulations  is  much  simpler  and  may  be  done  with  dependent  coordinates.  Therefore  the 
former  calls  for  special  purpose  implementations,  whereas  the  latter  is  more  suitable  for 


2.3  Refined  strategy  for  a  multi-index  variable  time  step  integration  method. 

A  refined  strategy  has  been  implemented  for  the  multi-index  variable  time  step  that  has  led  to  a 
better  behavior  and  control  of  the  time  step  size  during  the  integration  process  (publications  2 
and  5).  A  summary  of  the  method  and  results  are  reported  below,  (see  Appendix  for  more 

details). 

THE  MULTI  INDEX  AUTOSTEPPING  ALGORITHM 

As  reported  above  the  augmented  Lagrangian  index- 1  method  is  more  accurate  than  the  index-3 
for  small  time  steps,  and  the  index-3  method  is  more  efficient  than  the  index- 1  (with  same 
accuracy)  for  large  time  steps.  Therefore,  a  multi-index  fomiulation  would  take  advantage  of 
the  good  qualities  of  both  methods  for  different  time  step  sizes.  Index-3  is  the  nght  choice  for 
time  steps  between  10  ^  and  10 4  (where  it  is  accurate  and  efficient),  however,  as  we  decrease 
the  size  of  the  time  step,  starting  in  the  neighborhood  of  Af  =  10~4(with  no  pivoting)  or 
A;  =  10~5  (with  pivoting),  the  index- 1  method  becomes  more  accurate  and  robust. 

Also,  taking  into  account  the  correlation  existing  between  the  energy  and  the  local  error  an 
integral  of  motion  (system  total  energy  in  conservative  systems)  has  been  proposed  as  a 
measure  of  the  local  integration  error.  For  non-conservative  systems  and  considering  the  use  of 
fully  Cartesian  coordinates,  the  following  integral  of  motion  (energy  invariant)  is  established: 

n  =  r(0+v(0-J'qrQ^  =  constant 

where  T  is  the  kinetic  energy  and  V  the  potential  energy,  respectively.  This  equation  clearly 
yields  T+V  =  constant  for  conservative  systems. 

Using  this  energy  invariant  described  as  a  measure  of  the  local  integration  error,  a  variable 
time  step  size  strategy  has  been  proposed.  The  key  idea  is  to  modify  the  step  size  when  the 
change  in  the  energy  invariant  exceeds  an  allowable  value.  When  this  occurs  a  new  time  step 
size  is  calculated  in  order  to  achieve  the  desired  accuracy. 

NUMERICAL  RESULTS 

Simulations  are  performed  using  the  front  suspension  of  the  ILTIS  vehicle,  one  of  the  models 
included  in  the  library  of  benchmark  problems.  The  model  requires  a  total  of  23  variables, 
related  through  22  constraint  equations. 

The  simulation  is  performed  as  follows:  the  suspension  starts  moving  from  at  rest 
conditions  at  a  speed  of  5  m/s,  and  its  position  does  not  correspond  to  the  static  equilibrium. 
Thus,  it  freely  oscillates  until  the  static  equilibrium  position  is  reached.  After  three  seconds,  the 
suspension  goes  over  a  road  bump,  defined  by  a  cosine  function,  and  afterwards  it  freely 
oscillates  until  the  equilibrium  is  reached  again.  The  complete  analysis  lasts  for  6  seconds. 


general  purpose  simulators.  In  this  respect,  an  analogy  may  be  established  between  these 
methods  and  the  finite  element  formulation  for  structural  analysis. 

The  augmented  Lagrangian  method  (ALF)  features  a  leading  matrix  that  is  symmetric  and 
positive-definite.  This  makes  it  very  robust  in  simulations  with  singular  configurations  and 
redundant  constraints.  Also,  this  method  permits  the  differentiation  of  the  dynamic  equations 
to  get  the  tangent  matrix  in  a  Newton-Raphson  iteration,  which  allows  for  convergence  at 
larger  time  step  sizes  than  the  recursive  formulations.  Furthermore,  it  takes  a  serious 
advantage  of  parallel  computing.  However,  it  suffers  from  ill-conditioning  when  the  time 
step  size  becomes  small,  thus  having  convergence  difficulties  for  stiff  systems.  A  precision 
environment  of  32-digits  (words  of  128  bits)  could  be  the  cure  for  this  numerical  drawback. 

The  classical  Lagrangian  formulation  generates  an  almost  double-size  problem  than  its 
augmented  counterpart,  and  takes  less  profit  of  parallel  computing,  thus,  performing  notably 
worse.  In  addition,  its  leading  matrix  is  less  robust,  causing  the  failure  of  the  integration 
process  at  singular  configurations  and  in  the  presence  of  redundant  constraints.  However, 
this  method  is  not  so  sensitive  to  ill-conditioning  problems  for  small  time  step  sizes. 

The  fully-recursive  procedure  performs  well  for  stiff  and  non-stiff  systems.  Its  involved 
formulation  only  permits  a  fixed-point  iteration  scheme,  so  that  the  time  step  size  must  be 
kept  small  to  achieve  convergence.  Despite  of  this  fact,  it  shows  to  be  the  most  efficient 
method  among  the  three  compared,  but  does  not  obtain  any  advantage  from  parallel 
computing.  The  integration  also  fails  in  singular  configurations  and  in  the  presence  of 
redundant  constraints. 

According  to  the  presented  results,  real-time  simulation  of  large  multibody  systems  can  be 
achieved  with  medium-size  workstations,  opening  a  wide  field  for  development  of  hardware- 
in-the-loop  and  man-in-the-loop  facilities,  virtual  prototyping  tools  and  intelligent  vehicle 
control  systems. 


This  simulation  is  carried  out  under  two  different  conditions.  The  first  one  is  done  using  a 
index-3  approach  with  a  fixed  time  step  of  10~2  seconds.  The  second  simulation  is  performed  at 
variable  time  step  and  multi-index,  starting  with  index-3  and  a  step  size  of  10  3  seconds.  The 
first  simulation  took  8.41  seconds  of  CPU  time  and  the  second  took  36.07  seconds. 

Figure  6  shows  the  time  history  of  the  time  step  size.  It  may  be  clearly  seen  that,  for  the 
specified  tolerance  it  is  necessary  to  shorten  the  time  step  size  only  when  the  wheel  is  passing 
over  the  bump.  At  the  beginning  the  step  size  is  quickly  adjusted  to  speed  up  the  process  while " 
maintaining  the  local  error  between  the  allowable  limits. 


Simulation  time  (s) 

-  Index-3 
Index-1 

Figure  6.  Time  history  of  step  size. 


Figure  7  shows  the  total  system  energy  in  both  cases.  It  is  important  to  remark  that  only 
when  translating  the  origin  and  scaling  the  energy  represented  in  the  vertical  axis  one  can  notice 
its  change.  The  biggest  deviation  from  the  initial  value  (425,128  J  vs.  100  J)  is  about  0.02%  of 
the  total  system  energy  in  both  cases.  Also  important  it  is  the  fact  that  similar  precision  in  total 


system  energy  conservation  can  be  obtained  performing  the  simulation  with  an  index- 1 
approach  and  a  fixed  time  step  size  of  5.10'4  seconds.  In  that  case  the  simulation  takes  155 
seconds  (over  400%  more  than  in  the  variable  index-variable  step  case),  thus  corroborating  the 
benefits  of  using  the  proposed  muld-index  variable  time  step  procedure. 

As  a  result  of  this  example,  it  is  clear  that  a  multi-index- 1-3  approach  with  variable  time  step 
provides  an  excellent  solution,  since  the  more  efficient  index-3  method  is  used  during  a  large 
part  of  the  motion,  and  the  index- 1  only  at  that  part  where  the  time  step  decreases  below  the 
critical  limit  where  the  index-3  loses  precision  and  accuracy. 


Simulation  time(s) 

variable  time  step  size 
fixed  time  step  size 


Figure  7.  Total  system  energy. 


CONCLUSIONS 


•  The  system  energy  invariant  has  proved  to  be  a  good  measure  of  the  local  integration  error. 
This  feature  makes  unnecessary  the  use  of  traditional  (but  less  efficient)  error  criteria  used  in 
standard  numerical  integration  methods. 

•  The  variable  time  step  strategy  and  the  switch  from  one  method  to  the  other  do  not  cause 
problems  during  the  integration  process. 

•  The  use  of  the  proposed  technique  allows  a  substantial  speedup  gain  in  CPU  time  for  large 
scale  systems,  thus  helping  to  achieve  real-time  behavior. 

•  The  method  is  general  and  can  also  be  applied  to  solve  the  dynamics  of  flexible  multibodies. 
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Abstract.  Current  simulation  tools  for  multibody  dynamics  are  not  problem  dependent,  they  use  the 
same  modeling  process  to  all  cases  regardless  of  their  characteristics.  In  addition,  real-time  simulation 
of  small  multibody  systems  is  achievable  by  existing  simulation  tools,  however,  real-time  simulation 
of  large  and  complex  systems  is  not  possible  with  existing  methods.  This  is  a  challenge  that  needs 
to  be  addressed  before  further  advances  in  mechanical  simulation  with  hardware-in-the-loop  and 
man-in-the-loop,  as  well  as  virtual  prototyping  are  made  possible. 

This  paper  addresses  the  issue  of  how  the  modeling  process  -  dependent  versus  independent 
co-ordinates,  and  descriptor  form  versus  state-space  form  of  the  equations  of  motion  -  affects  the 
dynamic  simulation  of  multibody  systems  and  how  it  could  be  taken  into  account  to  define  the  concept 
of  intelligent  simulation.  With  this  new  concept  all  the  factors  involved  in  the  simulation  process  - 
modeling,  equations,  solution,  etc.  -  are  chosen  and  combined  depending  upon  the  characteristics 
of  the  system  to  be  simulated.  It  is  envisioned  that  this  concept  will  lead  to  faster  and  more  robust 
real-time  simulators. 

Key  words:  real-time-simulation,  computational  multibody  dynamics,  solution  algorithms. 


1.  Introduction  and  Background 

The  kinematics  and  dynamics  of  multibody  systems  constitute  an  important  part 
of  what  is  referred  to  as  CAD  and  MCAE.  The  advantage  of  computer  simulations 
performed  by  CAD  and  MCAE  tools,  is  that  they  allow  designers  to  predict  the 
kinematic  and  dynamic  behavior  of  all  types  of  multibody  systems  in  great  detail, 
during  all  the  design  stages;  from  the  first  design  concepts  to  the  final  prototypes. 
The  kinematic  and  dynamic  analysis  of  multibody  systems  are  processes  which 
are  most  appropriately  performed  in  real-time ,  allowing  even  the  analyst  to  act  as 
an  additional  element  in  the  simulation  process,  man-in-the-loop,  by  introducing 
external  or  control  forces  over  specific  degrees  of  freedom.  This  obviously  impos¬ 
es  constraints  on  the  computer  hardware  and  software.  Real-time  analysis  now 
requires  the  use  of,  at  least,  the  top  of  the  line  workstations,  and  is  not  yet  possible 
for  large  problems. 
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Figure  1  illustrates  all  the  factors  that  intervene  in  a  simulation  process.  Such 
factors  include:  modeling  process,  use  of  the  principle  of  mechanics,  formulation 
of  the  equations  of  motion,  numerical  integration  and  solution  of  the  resulting 
equations.  Intelligent  simulation  requires  the  right  choice  of  these  factors  which 
will  lead  to  the  best  possible  simulation  of  the  multibody  system  at  hand. 

The  first  issue  to  be  considered  in  the  simulation  is  that  of  modeling  the  system, 
that  is,  the  selection  of  a  set  of  parameters  or  co-ordinates  that  will  allow  to 
unequivocally  define  at  all  times  the  position,  velocity  and  acceleration  of  the 
system.  The  most  important  types  of  co-ordinates  currently  used  to  define  the 
motion  of  multibody  systems  are:  relative  co-ordinates,  reference  point  or  Cartesian 
co-ordinates,  and  natural  or  fully  Cartesian  co-ordinates.  A  thorough  comparison 
between  these  different  types  of  co-ordinates  and  their  influence  in  the  solution 
process  has  been  carried  out  by  Garcia  de  Jalon  and  Bayo  [1].  As  a  result  of 
that  study,  we  choose  in  this  paper  the  fully  Cartesian  co-ordinates  which  lead  to 
constant  inertia  terms  and  linear  Jacobian  of  the  kinematic  constraints,  and  are 
numerically  more  efficient  than  the  reference  point  co-ordinates  [1]. 
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Another  factor  in  the  modeling  process  is  the  choice  of  the  dy  namic  formulation 
obtained  from  the  application  of  the  principles  in  dynamics,  that  will  lead  o 
final  form  of  the  equations  of  motion.  Dynamic  principles  such  as  Lagrange  s  equa¬ 
tions,  Newton’s  laws,  canonical  equations  of  Hamilton,  virtual  power,  Hamilton  s 
principle,  Gibbs- Apell  equations,  etc.  constitute  the  basis  for  the  formula  ions  o 
multibody  dynamics.  The  choice  of  the  dynamic  formulation  will  influence  the 

subsequent  choice  of  numerical  integration  schemes. 

The  widely  used  method  of  Lagrange’s  multipliers  leads  to  a  represent  ion 
of  the  equations  of  motion  in  its  descriptor  form  constituting  a  set  of  index-3 
differential  algebraic  equations  (DAE).  The  addition  of  stabilization  techniques, 
such  as  the  method  of  Baumgarte  [2],  reduces  the  index  and  makes  the  solution 
tractable  by  means  of  standard  ODE  solvers,  however,  this  method  does  not  provide 
full  constraint  satisfaction,  leads  to  a  limited  control  of  accuracy,  and  in  add i  ion 
no  ways  are  available  for  choosing  the  values  of  the  coefficients  used  by  the 
method.  Recently,  an  augmented  Lagrangian  formulation  with  projections  [3]  was 
proposed  that  in  addition  to  transforming  the  set  of  equations  into  a  stabilize 
set  is  solvable  by  standard  ODE  methods  and  assure  Lyapunov  stability  of  the 
simulation  process  [4] .  They  also  have  the  advantages  of  being  robust  under  singular 
configurations,  topology  changes  and  with  redundant  constraints  and  provide  fu 

constraint  satisfaction.  ,  ,  rA1  A 

Other  formulations,  such  as  co-ordinate  partitioning  [5],  Kane  s  method  [6],  a 
virtual  power  with  projection  matrices  [7]  transform  the  equations  of  motion  to  a 
minimum  set  of  co-ordinates  or  state-space  form  that  is  directly  solvable  by  ODE 
methods.  State-space  representations  may  also  be  obtained  by  means  of  velocity 
transformations ,  initially  introduced  by  Jerkovsky  [8],  and  subsequently  extended 
by  other  authors  [9-13].  Recursive  methods,  also  called  topological  since they  t  e 
into  account  the  structure  of  the  mechanism,  such  as  those  developed  by  Walker 
and  Orin  [14]  or  Featherstone  [15,  16],  are  also  included  in  this  family.  State-space 
representations  are  more  suitable  for  ODE  integration  than  the  descriptor  counter¬ 
parts  at  the  expense  of  solving  the  velocity  and  position  problem  at  each  time  step. 
However,  they  do  not  handle  well  topology  changes  and  singular  configurations. 
In  addition,  they  cannot  support  stiff  integrators  and  consequently  may  present 

difficulties  when  the  system  has  built-in  numerical  stiffness. 

The  numerical  mathematics  community  has  sought  solutions  to  the  index  reduc¬ 
tion  problem  and  has  proposed  many  different  ways.  Recent  advances  have  been 
made  which  have  yielded  stabilized  index  reduction  methods  and  accurate  ways 
of  projecting  the  DAE  onto  the  underlying  ODE  for  more  stable  and  accurate 
solutions.  Key  developments  are  the  work  of  Brenan  et  al.  [17],  Gnepentrog  et  a  . 
[18],  Fiihrer  and  Leimkuhler  [19],  Lubich  et  al.  [20],  developer  of  MEXX,  Arnold 
[21]’  developer  of  HEDOPS,  Hairer  and  Wanner  [22],  developers  od  RADAU5, 

and  Petzold  [23],  developer  of  DASSL. 

In  this  paper  we  examine  and  compare  four  methods:  a  recently  propose 
augmented  Lagrangian  formulation  with  projections  [3]  in  index-1  and  index-3,  a 
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new  modified  state  space  formulation  based  on  projection  matrices,  that  is  capable 
of  handling  the  numerical  stiffness  coming  from  built  in  stiff  numerical  properties 
(springs,  dashpots,  etc.),  and  a  fully-recursive  formulation  proposed  by  Jimenez 
[24]  which  is  an  improvement  of  the  articulated  inertia  method  [15].  The  goal 
is  to  be  able  to  find  guidelines  as  to  what  modeling  methods  -  dependent  versus 
independent  co-ordinates,  descriptor  form  versus  state-space  form  of  the  equations 
of  motion,  global  versus  topological  methods  -  are  most  adequate  for  different  types 
of  multibody  systems.  For  this  reason,  the  comparison  is  made  through  the  analysis 
of  different  mechanical  systems  that  encompass:  open  and  close  chains,  systems 
that  undergo  singular  configurations,  modeling  with  redundant  constraints,  systems 
that  encounter  sudden  changes  and  high  frequency  thus  requiring  very  small  time 
steps  of  integration,  and  lastly  numerically  stiff  systems.  Comparisons  are  made  on 
the  basis  of  qualitative  as  well  as  quantitative  factors.  Conclusions  are  drawn  based 
on  the  performance  of  the  different  methods  that  serve  as  the  basis  to  establish  a 
strategy  for  intelligent  simulation. 

2.  Modeling  and  Dynamic  Analysis 

2.1.  Formulation  in  Fully  Cartesian  Co-ordinates 

Let  us  consider  a  multibody  system  whose  configuration  is  characterized  by  n 
fully  Cartesian  (or  natural)  co-ordinates  q  [1]  that  are  interrelated  through  the  m 
holonomic  kinematic  constraint  conditions 

3>(qJ,)  =  0.  (1) 

The  use  of  the  principle  of  virtual  power  directly  leads  to  the  equations  of  motion 

iqr(Mq-Q(q,q)  +  'tJ'A)  =0,  (2) 

which  for  a  general  multibody  system  leads  to: 

Mq  +  $[A  =  Q(q,q),  (3) 

where  M  is  the  mass  matrix,  Q  is  the  vector  that  contains  the  external  and  non¬ 
conservative  forces  as  well  as  the  velocity  dependent  inertia  forces,  3>q  is  the 
Jacobian  of  the  constraint  equations,  and  A  is  the  vector  that  contains  the  Lagrange 
multipliers.  Note  that  the  use  of  the  fully  Cartesian  co-ordinates  leads  to  a  constant 
mass  matrix  M  and  the  absence  of  velocity  dependent  inertia  forces  in  the  vector 
Q,  and  consequently  Equation  (3)  is  greatly  simplified. 

Equations  (1)  and  (3)  constitute  a  set  of  n  +  m  mixed  differential  algebraic 
equations  (DAE)  of  index-3  [17],  with  q  and  A  as  unknowns.  It  is  a  common 
practice  in  multibody  dynamics  to  double  differentiate  the  constraints  and  append 
the  resulting  equations  to  (3)  to  avoid  the  direct  integration  of  DAEs,  thus  yielding 
an  index- 1  formulation: 

‘  M  ]  f  q  1  f  Q 

<f>q  0  \  A  J  {  -<f>qq  -  $£ 
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These  equations  can  now  be  integrated  using  standard  numerical  integrators  with 
each  function  evaluation  performed  using  Equation  (4). 

2.2.  Index- 1  Augmented  Lagrangian  Formulation. with  Projections 

Recently,  a  method  was  presented  by  Bayo  and  Ledesma  [3]  that  uses  an  index- 1 
augmented  Lagrangian  method  with  mass-orthogonal  projections  of  the  positions 
and  velocities  to  their  constraint  manifolds.  Such  formulation  leads  to  the  following 
equations  of  motion: 

+  +  =Q(q,q),  (5) 

where  A*  are  the  Lagrange  multipliers.  Introducing  the  expression  $  —  $<jq  + 
+  <j>t  in  Equation  (5)  the  following  equation  is  obtained 

(M  +  <F'[a<L9)q  +  $^A*  =  Q(q,q)  -$^a(49q  +  $t)-  (6) 

It  is  important  to  note  that  there  is  a  substantial  difference  between  Equation  (6) 
and  the  Lagrange’s  multiplier  approach  represented  by  (4).  The  leading  matrix  of 
Equation  (4)  becomes  singular  in  singular  configurations,  changing  topologies  and 
in  the  presence  of  redundant  constraints.  However,  although  the  mass  matrixM  is  m 
general  positive  semi-definite,  the  leading  matrix  of  Equation  (6),  (M  +  <&q  a<&q), 
is  always  positive  definite,  which  means  that  it  can  always  be  factored,  even  in 
singular  positions,  changing  topologies  and/or  with  redundant  constraints. 

The  following  iteration  process  [3]  yields  the  unknown  multipliers  A 

A*+1  =  A*  +  i,  i  —  0, 1, 2, . . .  O) 

with  Aq  =  0  for  the  first  iteration.  Equation  (7)  physically  represents  the  introduc¬ 
tion  at  iteration  i  +  1  of  forces  that  tend  to  compensate  the  fact  that  the  addition  of 
all  the  constraint  terms  are  not  exactly  zero.  It  turns  out  that  with  the  augmented 
Lagrangian  formulation,  the  penalty  numbers  do  not  need  to  be  very  large  (thus 
leading  to  a  better  numerical  conditioning)  since  the  resulting  error  in  the  con¬ 
straint  equations  will  be  eliminated  by  the  Lagrange  terms  during  the  iteration 
procedure.  The  value  of  the  penalty  factor  a  does  not  affect  the  solution,  but  only 
the  convergence  rate.  Experience  shows  that  when  the  constraints  are  scaled  to 
unity,  penalty  factors  ranging  from  10s  to  I07  give  a  good  convergence  rate,  and 
only  1  to  2  iterations  are  required  to  converge  to  the  solution.  The  leading  matrix 

remains  constant  during  the  iteration  process. 

As  a  result  of  using  the  index- 1  formulation,  the  solution  of  Equation  (7)  yields 
a  set  of  accelerations  that  not  only  satisfies  dynamic  equilibrium  but  also  the  con¬ 
straint  conditions  $  -  0.  However,  it  is  expected  that  some  constraint  violations 
will  take  place  for  <E>  =  0  and  0.  A  projection  scheme  is  devised  for  q  and 
q  such  that  the  constraints  <£  and  $  will  be  zero  also.  These  projections  may  be 
done  as  follows: 
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(a)  Projection  in  q. 

During  the  time  integration  process  the  numerical  integration  scheme  yields  a  set  of 
co-ordinates  q*  that  does  not  completely  satisfy  the  constraint  conditions  $  =  0. 
In  order  .to  satisfy  the  constraints,  we  perform  a  mass-orthogonal  projection  of  the 
solution  to  the  constraint  manifold  and  obtain  a  new  set  of  positions  q  that  satisfies 
<f>  =  0.  This  can  be  achieved  by  the  solution  of  the  following  iterative  procedure 


[3]: 

(M  +  4>Ja$9)Aq(l+0  =  -7^ 

(8) 

with 

=  M(qW  -q*)  +  3^A(l) 

(9) 

and 

qfi+i)  =  q(0  +  Aq(l'+1V  A(j+I)  =  A(i)  +  a$(l+1), 

(10) 

where  the  superscript  indicates  the  iteration  number.  Equations  (8)  through  (10) 
can  be  used  iteratively  until  ||  Aq[|  <  e,  where  e  is  a  user  specified  tolerance. 

Note  the  important  feature  that  the  tangent  matrix  in  Equation  (8)  is  identical  to 
that  used  in  Equation  (6)  for  the  dynamic  analysis.  Furthermore,  since  the  solution 
q  is  close  to  the  initial  values  q*,  the  projection  problem  can  be  solved  using  a 
modified  Newton-Raphson  method  with  no  need  for  updating  the  tangent  matrix. 

(b)  Projection  in  q. 

Similarly  during  the  time  marching  integration  process  the  numerical  integration 
scheme  yields  a  set  of  velocities  q*  that  does  not  completely  satisfy  the  constraint 
conditions  4>  =  0.  Again,  we  perform  a  mass-orthogonal  projection  of  the  solution 
to  the  velocity  constraint  manifold  and  obtain  a  new  set  of  velocities  q  that  will 
satisfy  4>  =  0.  This  can  be  achieved  by  the  solution  of  the  following  equation  [3]: 


[M  -f  $[a$9]q  =  Mq*  -  +  cr),  (11) 

where  the  Lagrange  multipliers  of  the  projection  problem  are  updated  as: 

ff(i+i)=ff(i)+a§(i+|).  (12) 

A  more  efficient  implementation  of  the  projection  in  velocities  is  given  by  the 
following  recursive  set  of  equations: 

[M  +  #J’a#7]q(i+1)  =  Mq(i)-^J’a#t,  (13) 

and  the  iteration  is  started  by  setting 

q(°)=q*  and  <r(0)  =  (14) 


Since  the  leading  matrix  in  Equation  (11)  depends  only  on  the  generalized  co¬ 
ordinates,  it  needs  to  be  triangularized  only  once,  and  the  iteration  described 
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by  Equation  (13)  will  require  only  successive  forward  reductions  and  back-sub 
stitutions. 

2.3.  Index-3  Augmented  Lagrangian  Formulation  with  Projections 

In  [3],  an  index-3  augmented  Lagrangian  method  with  mass-orthogonal  projections 
of  the  velocities  and  accelerations  to  their  constraint  manifolds  is  presented.  Such 
a  formulation  leads  to  the  following  equations  of  motion. 

Mq  +  ^Q^  +  ^A*  =Q(q,q),  (15) 

where  A*  are  the  Lagrange  multipliers.  The  following  iteration  process  yields  the 
unknown  multipliers  A* 

A*+]  =  A*  +  ,  i  —  0, 1 , 2, . . .  (16) 

with  Aq  =  0  for  the  first  iteration.  Similar  to  the  case  of  the  index- 1  formulation  the 

value  of  the  penalty  factor  a  affects  the  convergence  rate.  Experience  shows  that 
when  the  constraints  are  scaled  to  unity,  penalty  factors  ranging  from  10  to  10 
give  a  good  convergence  rate  for  the  index-3  formulation,  and  only  1  to  2  iterations 

are  required  to  converge  to  the  solution. 

As  a  result  of  using  the  index-3  formulation,  the  solution  of  Equation  (15)  yields 
a  set  of  qn+,  that  not  only  satisfies  dynamic  equilibrium  but  also  the  constraint 
conditions  =  0.  However,  it  is  expected  that  some  constraint  violations  will  take 
place  for  4  =  0  and  =  0  .  A  projection  scheme  is  devised  for  q  and  q  such  that 
the  constraints  $  and  <£  will  be  zero  also.  The  projection  in  q  is  done  as  follows: 

(c)  Projections  in  q. 

Similar  to  the  velocity  analysis,  the  projection  of  the  generalized  accelerations 
onto  the  constraint  manifold  can  be  obtained  through  the  solution  of  the  following 
equation  [3]: 

(M  +  <&7q  a'i’Jq  =  Mq*  -  $q{a($Qq  +  #t)  +  «},  (^) 

where  the  associated  Lagrange  multipliers  are  updated  by  the  following  formula: 

*(*+1>  =  /c(i)  +  a<£(l+,).  (18) 

Numerical  implementation  of  the  augmented  Lagrangian  formulation  for  the  pro¬ 
jection  of  accelerations  is  given  by  the  recursive  set  of  equations 

(M  +  $Ja$q)q(i+1)  =-  Mq(t)  -  O9) 

and  the  recursion  is  started  by  setting 

q(°)  =  q*  and  k, (°)  =  a$(0^. 


(20) 
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Equation  (19)  constitutes  a  system’ of  linear  algebraic -equations  in  q  that  are 
solved  iteratively  until  convergence  in  the  generalized  accelerations  is  achieved. 
Since  the  leading  matrix  in  (19)  is  the  same  as  that  of  Equation  (13)  used  for  the 
velocity  projection,  the  acceleration  analysis  will  require  only  successive  forward 
reductions  and  back-substitutions. 

2.4.  Modified  State-Space  Formulation 

A  state-space  representation  of  the  equations  of  motion  can  be  obtained  using  the 
basis  of  the  null-space  of  the  Jacobian.  Such  basis  is  defined  by  the  matrix  R  which 
will  satisfy  the  following  condition  [1]: 

<^R  -  0.  (21> 

Accordingly,  the  dependent  velocities  and  accelerations  can  be  expressed  in  terms 
of  the  independent  ones  as  follows: 

q  =  Rz  and  q  =  Rz  +  Rz  (22) 

Introducing  Equation  (22)  into  Equation  (3)  the  following  equation  is  easily 
obtained: 

RtMRz  =  RTQ(q.  q)  -  RTMRz,  (23) 

which  constitutes  the  equations  of  motion  in  independent  co-ordinates.  Equation 
(23)  provides  the  function  evaluation  for  the  numerical  integration  scheme.  This 
method  requires  the  solution  of  the  position  and  velocity  problems  [1]  at  each  time 
step  to  obtain  the  matrix  R,  q  and  q,  so  that  Equation  (23)  may  be  solved.  The 
solution  of  the  position  and  velocity  problems,  that  are  costly  computationally, 
have  the  same  meaning  as  that. of  projecting  the  solution  to  obtain  the  dependent 
positions  and  velocities. 

Implicit  integrators  make  successive  calls  to  (23)  as  a  fixed  point  iteration. 
Such  iteration  is  either  generally  inefficient  or  precludes  convergence  for  stiff 
systems.  Due  to  this  reason  we  propose  below  a  modification  of  this  method 
based  on  the  direct  combination  of  Equation  (23)  with  the  difference  equations  of 
the  numerical  integration  scheme.  Such  combination  permits  a  modified  Newton- 
Raphson  iteration  that  allows  this  method  to  be  used  in  stiff  problems,  and  is 
described  in  detail  below. 

2.5.  Fully-Recursive  Formulation 

A  topological,  fully-recursive  algorithm  of  order  O(N)  has  been  presented  by 
Jimenez  [24],  This  algorithm  uses  a  set  of  relative  co-ordinates  z  to  represent  the 
state  of  the  mechanism  and  the  Cartesian  velocities 

z‘  =  {4 


(24) 
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and  their  derivatives  Zj  to  formulate  the  equations  of  motion.  The  vector  roi 
represents  the  velocity  of  the  point  of  body  i  that  is  currently  situated  at  the  origin 
of  the  inertial  reference  frame,  while  u>l  is  the  angular  velocity  of  body  i.  In  the 
case  of  open  chain  mechanisms,  the  relative  co-ordinates  z  are  independent,  equal 
in  number  to  the  degrees  of  freedom.  If  there  are  close  chains,  they  are  opened 
either  at  a  pair  or  body,  and  the  corresponding  constraints  are  established. 

Once  the  co-ordinates  are  defined,  the  first  step  consists  of  a  kinematic  trans¬ 
mission  starting  from  the  base  towards  the  end  of  the  chain.  The  corresponding 
recursive  expressions  in  velocities  and  accelerations  are  given  by 


Zl  =  Z%-\  +  bjZ;, 


Z  i  =  Z,_1  +  biit  +  d1; 


(25) 


(26) 


where  bz  is  a  6  x  1  vector  containing  the  value  of  the  Cartesian  velocities  Z;  when 
all  the  velocities  z  are  made  zero  except  Z;  =  1,  and  d;  is  also  a  6  x  1  vector 
containing  the  difference  of  the  Cartesian  accelerations  Z*  -  Zj_i,  when  all  the 
accelerations  z  are  made  zero. 

The  second  step  is  a  condensation  of  the  forces  and  inertial  terms  starting 
from  the  end  of  the  chain  and  progressing  towards  the  base.  In  cases  with  closed 
chains,  the  reaction  forces  coming  from  points  where  chains  have  been  opened 
are  introduced  through  a  penalty  formulation  as  functions  of  the  violation  of  the 
constraints  and  their  derivatives.  These  reaction  forces  are  condensed  together  with 
the  rest  of  the  forces.  The  recursive  expressions  for  this  second  step  are 


Mi-i  =  Mi_i  +  KjMj, 

Qi_i  =  Qi_i  +  Ki(Qt  -  Mldi), 


(27) 

(28) 


where 


Ki  =  U 


Mtbjbj 

b-Mjbj 


(29) 


Mi  and  Qz  are,  respectively,  the  mass  matrix  and  the  vector  of  generalized  forces 
of  body  i. 

The  third  and  last  step  is  the  calculation  of  accelerations  z  from  the  base  to  the 
end  of  the  chain.  For  this  purpose,  the  following  expressions  are  provided 


z,  - 


Z,;  = 


bj(Qi-M.di) 
b£,  M  i  b  i 

b£[QI-Mi(dl  +  Zi_1)] 

b-Mibj 


(30) 


(31) 
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Therefore,  it  may  be  seen  that  the  formulation  is  fully  recursive  and  there  is  no  need 
of  solving  a  system  of  equations.  This  fact  suggests  a  good  speed  competitiveness 
of  this  method  as  the  size  of  the  problem  increases.  In  those  cases,  non-recursive 
formulations  are  forced  to  cope  with  large  sets  of  equations.  On  the  other  hand  the 
complexity  of  this  formulation  and  its  implementation  is  far  superior  to  those  seen 
above,  especially  the  index- 1  and  index-3  methods. 

3.  Numerical  Implementation  of  the  Different  Methods 

The  difference  equations  of  the  numerical  integration  scheme  may  be  expressed  in 
general  terms  as: 

<Wi  =  aqn+l  -  qn  and  qn+l  =  &qn+1  _  <jn>  (32) 

where  a  and  b  are  constants  that  depend  on  the  numerical  integration  scheme  and  the 
time  step,  qn  and  qn  are  known  quantities  that  depend  on  the  positions,  velocities 
and  accelerations  at  step  n  and/or  previous  steps.  Note  that  if  the  method  is  explicit 
a  and  b  will  be  equal  to  zero. 


3. 1 .  Index-1  Augmented  Lagrangian  Formulation 
Substituting  Equation  (32)  in  the  equation  of  motion  (5)  for  index- 1,  yields: 
(M  +  *?CL*q)bqn+x  +  +  *rQ(49  q  +  *t) 


-  Q(q,  q)  -  (M  +  $>^a$>q)q  =  o,  (33) 

which  constitutes  a  set  of  n  equations  of  the  form  f  =  0  to  be  solved  for  the 
unknown  position  qn+I.  We  can  apply  a  modified  Newton-Raphson  method  and 
evaluate  the  solution  by  means  of  the  iterative  process: 


■<9fl  W 

dq. 


AqF+0  =  _f(q(0)j 


(34) 


where  the  function  f  is  defined  as 


f(ln+l)  =  Mq„+I  +  +  A*)  -  Qn+I, 

and  the  tangent  matrix  is  approximated  by  the  following  matrix: 

Q^  =  HM  +  ^^q)-Qq-aQq, 


(35) 


(36) 


where  the  last  two  terms  of  (36)  represent  the  contribution  to  the  tangent  matrix 
elast,c  forces  ('-e-  sPrings)  and  those  dependent  on  the  velocity 
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3.2.  Index-3  Augmented  Lagrangian  Formulation 

The  substitution  of  Equation  (32)  into  Equation  (15)  for  index-3  yields 

£>Mqn+j  +  <F^(ck3?  +  A*)  =  Qn+i  +  @7) 

which  again  constitutes  a  set  of  nonlinear  algebraic  equations  with  qn+1  as 
unknowns.  We  use  a  modified  Newton-Raphson  iteration  and  evaluate  the  solution 
by  means  of  the  iterative  process  outlined  above  with  the  function  f  defined  as 

f(q„+i)  =  Mqn+1  +  +  A*)  -  Qn+i  (38) 

and  the  tangent  matrix  approximated  by 


-  =  6M  +  ~Qq~  aQq- 

oq 


(39) 


3.3.  New  State-Space  Implementation  for  Stiff  and  Non-Stiff  Systems 

One  of  the  major  drawbacks  of  the  state-space  representation  given  by  Equation  (23) 
is  that  this  equation  is  customarily  used  as  a  fixed  point  iteration  for  function 
evaluation.  It  is  well  known  in  numerical  analysis  that  such  iteration  is  not  efficient 
and  even  fails  when  the  equations  of  motion  are  stiff.  Consequently,  in  what  follows, 
we  propose  a  combination  of  the  equations  of  motion  (23)  and  the  finite  difference 
equations  of  the  numerical  integration  scheme  (32),  such  that  the  integration  of  both 
stiff  and  non-stiff  problems  can  be  made  tractable.  Since  we  are  now  integrating 
using  independent  co-ordinates,  Equation  (32)  may  be  rewritten  as: 

zn+i  =  azn+i  —  zn  and  zn+i  =  bzn+\  —  zn.  .  (40) 

The  substitution  of  Equation  (40)  in  the  equation  of  motion  for  state-space  (23), 
yields: 

RTMR6zn+l  -  RTQ(q,q)  -  RTMRazn+1  -  RTMRi  -  RTMRi  =  0,(41) 

which  constitutes  a  set  of  n  equations  of  the  form  f  =  0  to  be  solved  for  the 
unknowns  zn+i.  Again,  we  can  apply  a  modified  Newton-Raphson  method  and 
evaluate  the  solution  by  means  of  the  iterative  process. 


dfi(l) 

dz 


Az<i+.) 


where  the  function  f  is  defined  as: 

f(zn+i)  =  RrMRzn+1  -  RTQ7l+  i  +-RTMRzn+1  =  0, 
and  the  tangent  matrix  is  approximated  by  the  following  matrix: 


(42) 


(43) 


270 


J.  CUADRADO  ET  AL. 


Table  I.  Comparative  results  for  the  double  pendulum. 


A  t  = 

1(T2 

At  = 

10~3 

At  = 

10~4 

At  = 

10~5 

Error 

Time 

Error 

Time 

Error 

Time 

Error 

Time 

index-1 

No  convergent. 

1.1E-3 

2.8 

1 . 1 E-5 

21.4 

2.5E-7 

214 

Index-3 

2.1E-1 

0.4 

4.6E-3 

1.9 

2.0E-5 

17.2 

Wrong  results 

State-space 

2.9E-1 

1.8 

4.8E-3 

10.5 

4.3E-5 

96.6 

7.6E-7 

862 

Fully-rec 

2.7E-1 

2.0 

2.6E-3 

11.3 

3.4E-5 

100.3 

3.4E-7 

1003 

where  the  last  two  terms  of  (44)  represent  the  contribution  to  the  quasi-tangent 
matrix  coming  from  the  elastic  forces  (i.e.  springs)  and  forces  dependent  on  the 
velocity  (i.e.  dashpots).  The  contribution  of  elastic  and  damping  terms  to  the 
tangent  matrix  makes  this  method  more  suitable  and  efficient  for  stiff  systems  than 
the  classical  fixed  point  iteration  which  only  includes  the  first  term  in  the  r.h.s.  of 
(44). 

3.4.  Fully-Recursive  Formulation 

As  it  may  be  easily  understood  after  the  description  of  this  formulation  given  in 
the  previous  section,  the  own  nature  of  the  fully  recursive  method  makes  difficult 
to  introduce  the  integrator  scheme  (40)  into  the  dynamic  equations,  since  these 
equations  are  not  explicitly  formulated,  but  developed  in  three  different  recursive 
steps.  Therefore,  in  this  case  the  integration  of  the  equations  of  motion  must  be 
performed  using  a  fixed  point  iteration.  The  corresponding  function  evaluations 
(acceleration  analysis)  are  provided  by  the  described  recursive  formulation,  and 
the  integration  scheme  (40)  changes  its  form  to 

zn+)  =  |  zn+i  -  iTl  and  in+\  =  -  zn+i  -  zn.  (45) 

b  cl 

This  fact  constitutes  a  drawback  when  comparing  this  method  to  the  other  formu¬ 
lations  described  above  which  apply  a  Newton-Raphson  scheme.  However,  as  said 
before,  this  method  does  not  require  the  solution  of  any  set  of  linear  equations.  We 
will  show  below  which  effects  dominate  in  different  cases. 

4.  Numerical  Results 

First  Example.  The  simulation  of  a  double  pendulum  that  moves  from  rest  in  the 
horizontal  position  under  gravity  effects  is  carried  out.  Each  link  has  a  distributed 
unit  mass  and  a  unit  length .  We  use  natural  co-ordinates  for  the  modeling  process  [  1  ] 
and  for  the  integration  we  use  the  trapezoidal  rule  which  is  second  order,  implicit, 
A-stable  method  and  also  energy  preserving  in  the  linear  regime.  The  total  time  of 
simulation  is  10  seconds  and  the  penalty  factor  is  107. 

Table  I  shows  the  maximum  error  in  the  energy  norm  and  the  CPU  time  in 
seconds  taken  by  each  of  the  four  methods  for  different  time  steps  of  integration. 
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Figure  2.  Assembly  of  two  four-bar  linkages. 


It  is  important  to  emphasize  that  the  index-3  formulation  does  not  converge  for  the 
time  steps  1CT4  and  10-5  unless  the  penalty  factor  is  increased  to  109  and  1011, 
respectively.  Even  in  the  last  case  the  results  obtained  are  wrong  due  to  the  ill- 
conditioning  caused  by  the  small  time  steps  in  the  index-3  formulation  (the  reader 
is  referred  to  [17]  for  an  explanation  of  this  effect).  Note,  however,  that  no  pivoting 
or  filtering  is  done  in  the  solution  process  to  avoid  or  diminish  this  ill-conditioning. 

The  index- 1  method  is  more  robust  than  the  index-3.  It  does  not  require  to 
modify  the  penalty  factor  and  increases  its  accuracy  as  the  time  step  is  decreased, 
as  clearly  seen  in  the  table.  The  index-3  method  is  more  efficient  than  the  index- 1 
for  large  time  steps  of  integration. 

The  state-space  and  the  fully  recursive  methods  work  for  all  the  time  steps 
considered.  Therefore,  they  show  a  robust  behavior  although  they  run  notably 
slower  than  its  two  competitors.  In  this  example,  the  slowness  of  the  fourth  method 
is  justified  due  to  the  reduced  number  of  variables  used  by  the  other  methods 
(just  4).  Moreover,  the  fully-recursive  formulation  has  been  developed  for  three- 
dimensional  cases  and  6x6  matrices  and  6  x  1  vectors  are  used.  Undoubtedly,  this 
formulation  could  be  simplified  for  planar  cases,  thus  getting  smaller  calculation 
times. 

Second  Example.  The  simulation  of  a  one  degree-of-freedom  assembly  of  two 
four-bar  linkages  has  been  proposed  in  [25]  as  an  example  to  test  the  performance 
in  cases  where  the  mechanism  undergoes  singular  configurations.  When  the  mech¬ 
anism  reaches  a  horizontal  position,  the  number  of  degrees  of  freedom  increases 
instantaneously  from  1  to  3.  To  define  the  position  of  the  system,  we  use  the  six 
position  variables  (aq,  y i,  xj,  yi,  £3,  2/3),  as  may  be  seen  in  Figure  2.  All  the 
links  are  of  length  l  =  1  m  and  have  a  uniformly  distributed  mass  m  =  1  kg. 
The  gravity  force  acts  in  the  negative  Y  direction,  with  a  value  g  —  9.81  m/s2.  At 
i  =  0  the  initial  velocity  is  x  1  —  1.  We  integrate  the  motion  for  10  seconds  with 
107  the  value  of  the  penalty  parameter.  Figure  3  shows  different  time-histories  of 
the  simulation  and  shows  that  the  mechanism  goes  through  the  singular  position 
(|xi  |  —  1)  10  times. 
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Figure  3.  Time-history  of  the  two  four-bar  linkages. 


The  results  obtained  are  illustrated  in  Table  II.  The  same  trend  observed  in  the 
previous  case  occurs  again  in  this  one.  The  indcx-1  method  is  more  robust  in  the 
sense  that  it  increases  its  accuracy  as  A t  decreases  although  it  does  not  converge  for 
big  time  steps.  The  index-3  is  more  efficient  and  converges  for  bigger  time  steps 
but  shows  moderate  errors  for  10-3  and  fails  for  time  steps  equal  to  or  smaller 
than  10~4.  As  expected  [25],  neither  the  state-space  method  nor  the  fully-recursive 
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Table  II.  Comparative  results  for  the  two  four-bar  linkages. 


At  =  5  *  1(T2 

At  = 

10“2 

II 

HO 

< 

10~3 

At  = 

10~4 

Error  Ti  me 

Error 

Time 

Error 

Time 

Error 

Time 

Index- 1 

No  convergent. 

5.6E-3 

1.0 

3.0E-4 

7.2  . 

2.8E-7 

63.6 

Index-3 

7.6E-2  0.3 

5.0E-3 

0.6 

7.7E-2 

4.0 

Wrong  results 

State-space 

Does  not  work 

Does  not  work 

Does  not  work 

Does  not  work 

Fully-rec 

Does  not  work 

Does  not  work 

Does  not  work 

Does  not  work 

(x4.y4) 


Figure  4.  Andrew’s  mechanism. 

formulation  work  in  these  kind  of  simulations  because  they  either  fail  or  lead  to 
divergent  results  when  the  singular  position  is  reached. 

Third  Example.  The  Andrew’s  squeezing  mechanism,  described  in  detail  in  [26], 
has  been  proposed  as  an  example  to  test  numerical  algorithms  in  multibody  dynam¬ 
ics.  It  requires  very  small  time  steps  because  of  the  small  time  scale  of  the  problem. 
We  perform  the  simulation  using  natural  co-ordinates  and  redundant  constraints 
to  introduce  the  angle  of  the  crank  (3  as  an  additional  variable  (see  Figure  4).  The 
penalty  factor  is  a  —  1 07,  and  again  the  trapezoidal  rule  is  used  for  the  numer¬ 
ical  integration  for  a  total  time  of  0.03  seconds.  Sparse  solvers  are  used  in  the 
non-recursive  methods  to  solve  the  system  of  equations.  The  simulation  yields  the 
values  of  the  angles  shown  in  Figure  5,  which  are  calculated  in  a  post-process  as 
functions  of  the  natural  co-ordinates  used  to  model  the  mechanism.  The  angles  in 
Figure  4  and  the  plot  numbers  in  Figure  5  are  related  as  follows:  1  —  (3,  2  —  9,  3  -  7, 
4  —  8, 5  —  <f>,  6  —  Q.  The  analysis  is  performed  using  the  four  proposed  methods  and 
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Figure  5.  Angles  of  the  Andrew’s  mechanism. 

Table  III.  Comparative  results  for  the  Andrew’s  mechanism. 


At  = 

1(T4 

At  - 

10"5 

At  = 

10~6 

Error 

Time 

Error 

Time 

Error 

Time 

Index- 1 

1 .4E-2 

0.6 

2.5E-2 

5.2 

2.0E-4 

54 

Index-3 

Wrong  results 

No  convergent. 

No  convergent. 

State-space 

4.3E-3 

1.3 

9.3E-5 

10.2 

6.2E-5 

101.9 

Fully-rec 

3.0E-3 

2.8 

9.3E-5 

16.2 

3.2E-5 

145.6 

the  results  are  depicted  in  Table  III.  It  may  be  seen  that  since  the  time  step  required 
is  very  small,  the  index-3  method  either  fails  or  yields  wrong  results  (again  no 
pivoting  or  filtering  is  performed).  On  the  other  hand,  the  index- 1  method  yields 
good  results  and  shows  better  accuracy  as  the  time  step  decreases.  The  state-space 
method  shows  this  behavior  also,  giving  even  a  slightly  better  energy  conservation 
than  index- 1,  but  running  at  half  of  the  speed.  The  fully-recursive  formulation 
presents  also  a  robust  behavior,  although  it  runs  the  slowest  of  all  because  of  the 
number  of  close  kinematic  loops  that  the  problem  has.  The  system  of  equations  for 
the  index- 1  and  index-3  methods  are  solved  using  solvers  that  take  advantage  of 
the  sparse  nature  of  their  leading  matrices. 

Fourth  Example.  The  4  x  4  litis  vehicle  [27]  has  been  proposed  as  a  benchmark 
problem  by  the  European  automobile  industry  to  check  multibody  dynamic  codes. 
We  perform  two  different  simulations  using  the  front  suspension  of  the  vehicle  (see 
Figure  6).  The  model  requires  a  total  of  23  variables,  related  through  22  constraint 
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equations  with  only  one  degree  of  freedom.  The  characteristics  of  the  suspension 
are: 

—  A  shock  absorber,  which  provides  a  damping  force  given  by  the  following 
nonlinear  expression: 

Fd  =  9945.627v  +  33955.72v2  -  59832.25v3  -  39565 l.Ov4  [N] 
for  —0.2  <  v  <  0.21  m/s 

Fd  =  -416.42  +  1844. 3v  [N]  forv<-0.2m/s 

Fd  =  1919.1638+  1634.727v  [NJ  for  v  <  0.21  m/s 

and  a  elastic  force  due  to  an  external  polymer,  given  by 

Fs  =  -4.0092  *  106  +  2.8397  *  107x 

-  6.7061  *  10V  +  5.2796  *  10V  [N], 

where  the  distance  x  is  defined  in  meters. 

—  A  tire,  modeled  by  means  of  a  linear  vertical  spring  of  stiffness  460,000  N/m. 

—  A  leaf  spring,  modeled  as  a  linear  spring  of  35,900  N/m  for  deformations 
smaller  than  14.5  cm.  When  such  deformation  is  exceeded  the  suspension  hits 
a  second  spring  of  value  107  N/m,  thus  adding  a  much  higher  stiffness  to  the 
suspension. 

The  first  simulation  is  performed  as  follows:  the  suspension  starts  moving  from 
at  rest  conditions  at  a  speed  of  5  m/s,  and  its  position  does  not  correspond  to  the 
static  equilibrium.  Thus,  it  freely  oscillates  until  the  static  equilibrium  position  is 
reached.  After  three  seconds,  the  suspension  goes  over  a  road  bump,  defined  as  a 
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Table  IV.  Comparative  results  for  the  litis  front  suspension. 


At  = 

1(T2 

At  = 

10'3 

II 

< 

10"4 

At  = 

10“5 

Error 

Time 

Error 

Time 

Error 

Time 

Error 

Time 

Index- 1 

No  convergent. 

1.2E-3 

101.2 

6.0E-5 

464.9 

3.0E-7 

4791 

Index-3 

1.6E-1 

6.0 

2.1E-3 

38.3 

7.5E-2 

377.2 

Wrong  results 

State-space 

1.6E-1 

16.7 

7.0E-4 

137.0 

4.0E-5 

1091.0 

4.0E-7 

10435 

Fully-rec 

6.6E-2 

6.6 

4.2E-4 

38.6 

6.0E-6 

392.0 

2.0E-7 

4079 

Figure  7.  Vertical  position  of  the  wheel  center.  Non-stiff  case. 


cosine  function,  and  afterwards  oscillates  until  the  equilibrium  is  reached  again. 
The  stiff  spring  of  value  107  N/m  is  not  activated  in  this  simulation,  therefore,  it 
cannot  be  considered  as  a  stiff  problem.  The  complete  analysis  lasts  for  6  seconds. 
The  analysis  is  carried  out  using  the  four  methods  and  the  performance  results  are 
depicted  in  Table  IV.  Figure  7  illustrates  the  time  history  of  the  center  of  the  wheel. 

Table  IV  shows  that  the  index-3  method  does  not  converge  for  time  steps  10~4 
and  10-5  unless  the  penalty  factor  is  increased  to  1011  and  1013,  respectively.  For 
a  step  of  10"5  the  results  obtained  are  wrong  because  of  the  ill-conditioning.  The 
index- 1  method  is  again  much  more  robust,  not  requiring  to  modify  the  penalty 
factor,  whose  value  is  set  to  107  for  all  the  time  steps.  Its  accuracy  increases  as 
the  time  step  is  decreased,  but  the  method  does  not  converge  for  A t  —  10~2  (no 
pivoting  is  used).  The  state-space  method  also  shows  a  robust  behavior  although, 
in  this  case,  it  runs  notably  slower  than  the  other  methods. 

Table  IV  also  shows  that  the  fully-recursive  formulation  becomes  competitive  in 
the  simulation  of  this  non-stiff  large  system.  The  reason  for  this.behavior  is  that  in 
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this  example  the  first  non-recursive  methods  use  23  variables,  notably  more  than  in 
the  previous  examples,  and  require  the  triangulanzation  of  matrices  of  size  23  *  23 
(sparse  matrix  algebra  has  been  used  to  solve  these  equations).  Conversely,  the 
fully-recursive  method  does  not  require  the  solution  of  such  system  of  equations 
in  each  function  evaluation.  It  is  important  to  note  how  the  index-3  method  and  the 
fully-recursive  method  take  very  similar  times  of  integration,  thus  corroborating 
the  equivalence  between  exploiting  sparsity  and  the  articulated  body  approach  for 
systems  with  a  large  number  of  bodies. 

The  second  simulation  is  performed  activating  the  second  spring  so  that  when 
the  suspension  hits  the  bump  and  the  leaf-spring  exceeds  a  deformation  of  14.5  cm 
the  spring  of  value  10  N/m  becomes  active,  thus  introducing  additional  stiffness  to 
the  system.  The  simulation  is  performed  using  the  index- 1  and  the  fully  recursive 
formulations  at  different  time  steps  of  integration.  Figure  8  shows  the  time  history 
o  the  center  of  the  wheel  for  this  case.  It  is  important  to  point  out  that  the  fully- 
recursive  method  requires  a  much  smaller  time  step,  5.0*  10~6  seconds  compared  to 
the  index- 1  method  which  successfully  performs  the  simulation  starting  at  5.0*  10~5 
seconds.  It  becomes  obvious  that  for  stiffer  systems  such  as  those  resulting  from 
■argcr  springs,  visco-elastic  bushings,  flexible  multibodies,  active  control  systems 
etc.,  the  performance  of  the  fully  recursive  method  will  further  degrade. 

As  a  result  of  this  example,  it  is  clear  that  a  multi-index- 1-3  approach  with 
variable  time  step  would  provide  an  ideal  solution,  since  the  more  efficient  index-3 
method  would  be  used  during  a  large  part  of  the  motion  and  the  index- 1  only  at 
that  part  where  the  time  step  decreases  below  10~4. 
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5.  Conclusions 


We  have  addressed  in  this  paper  the  issue  of  how  the  modeling  aspects  of  depen¬ 
dent  versus  independent  co-ordinates,  descriptor  form  versus  state-space  form  of 
the  equations  of  motion,  and  global  versus  topological  methods  affect  the  dynamic 
simulation  of  multibody  systems.  More  concretely,  we  have  examined  and  com¬ 
pared  the  use  of  the  top  of  the  line  methods  that  include:  an  augmented  Lagrangian 
formulation  with  projections  in  index- 1  and  index-3 ,  a  new  modified  state-space 
formulation  based  on  projection  matrices,  that  is  capable  of  handling  the  numerical 
stiffness  coming  from  built  in  stiff  numerical  properties  (springs,  dashpots,  etc.), 
and  a  fully -re  cursive  formulation  whose  calculation  effort  grows  linearly  with  the 
size  of  the  problem.  The  aim  is  to  assess  how  these  factors  may  be  taken  into 
account  to  establish  a  new  concept  in  intelligent  simulation. 

These  are  the  major  conclusions: 

•  The  index-3  formulation  with  projections  is  the  most  efficient  of  all  the  meth¬ 
ods  tested,  but  it  leads  to  serious  ill-conditioning  at  small  time  steps,  starting  at 
about  10~5,  when  the  penalty  parameter  becomes  problem  dependent.  (Note 
that  neither  pivoting  nor  filtering  techniques  have  been  used  to  avoid  or  dimin¬ 
ish  ill-conditioning.) 

•  The  index- 1  formulation  with  projections  is  the  most  robust  of  all  the  methods, 
but  it  cannot  converge  for  large  time  steps  where  the  index-3  does.  Contrary 
to  index-3  it  becomes  more  accurate  as  the  time  step  decreases  and  it  is  not 
affected  by  ill-conditioning.  It  is  also  insensitive  to  the  penalty  parameter 
which  remained  constant  in  all  the  simulations. 

•  The  new  state-space  method  works  well  with  stiffness,  shows  robust  behavior 
and  good  energy  preservation,  as  well  as  good  convergence  at  all  time  steps. 
However,  it  runs  considerably  slower  than  the  index- 1  and  index-3  methods  in 


all  the  cases  tested. 

•  The  fully-recursive  (articulated  inertia)  formulation  behaves  in  a  similar  way 
to  that  of  the  previous  one,  except  in  what  refers  to  speed  and  stiffness.  This 
method  is  much  faster  than  the  state-space  but  not  suitable  for  stiff  problems 
(such  as  those  resulting  from  large  springs,  visco-elastic  bushings,  flexible 
multibodies,  active  control  systems,  etc.).  Also,  it  turns  out  that  for  small 


systems  its  performance  is  inferior  to  the  index-3  and  index- 1  methods,  but 
improves  as  the  size  of  the  problem  increases.  For  large  non-stiff  problems  the 
method  becomes  competitive. 

The  following  tables  summarize  the  findings  of  this  paper.  Each  method  is 
given  a  grade  depending  on  its  performance  under  different  situations.  The 
letter  “A”  means  outstanding,  “B”  good,  “C”  fair,  “D”  poor  and  “F”  signifies 


a  failure. 
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Method 

Easiness 
of  implemen¬ 
tation 

Free 

from  graph 
methods 

Speed 

Accuracy 

Small 
to  medium 
problems 

Medium 
to  large 
problems 

At  <  nr5 

At  >  lO-5  At  >  10~2 
At  <  icr2 

Index- 1 

A 

A 

B 

B 

A 

A  F 

lndex-3 

A 

A 

A 

A 

F 

A  A 

Space-state 

B 

A 

C 

C 

A 

A  A 

Recursive 

C 

D 

C 

A 

A 

A  A 

Method 

Performance 

Changing 

topologies 

Singular 

configurations 

Inequality 

constraints 

Numerical 

stiffness 

Redundant 

constraints 

Index- 1 

A 

A 

A 

A 

A 

Index-3 

A 

A 

A 

A 

A 

Space-state 

C 

F 

C 

A 

C 

Recursive 

D 

F 

D 

D 

D 

•  It  remains  to  address  in  future  work  the  real-time  behavior  of  these  methods  on 
large  scale  industrial  problems,  such  as  the  complete  model  of  an  automobile. 
More  concretely,  a  multi-index- l/index-3  formulation  could  lead  to  an  ideal 
situation  since  it  encompasses  the  benefits  of  both  methods  [28],  Also  sparse- 
parallel  solvers  have  a  lot  of  potential  for  the  index- l/index-3  approach  since  a 
great  deal  of  operations  may  be  performed  element-by-element.  On  the  other 
hand,  the  fully  recursive  formulation  becomes  competitive  on  large  non-stiff 
problems,  although  to  be  considered  as  a  general  tool  it  will  require  graph 
methods  for  identifying  close-chain  loops  and  spanning  trees.  Therefore,  the 
recursive  method  is  a  good  candidate  for  special  purpose  simulators  of  large 
non-stiff  multibody  systems. 
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SUMMARY 

This  paper  presents  a  multi-index  variable  time  step  method  for  the  integration  of  the  equations  of  motion  of 
constrained  multibody  systems  in  descriptor  form.  The  basis  of  the  method  is  the  augmented  Ugrang.an  formulation 
with  projections  in  index-3  and  index- 1.  The  method  takes  advantage  of  the  better  performance  of  the  index-3 
formulation  for  large  time  steps  and  of  the  stability  of  the  index-1  for  low  time  steps,  and  automatically  switches  from 
one  method  to  the  other  depending  on  the  required  accuracy  and  values  of  the  time  step.  Various  numerical  problems  that 
arise  during' the  simulation  process  are  described.  The  paper  also  describes  ways  to  circumvent  problems. 

The  variable  time  stepping  is  accomplished  through  the  use  of  an  integral  of  motion,  which  in  the  case  of  conservative 
systems  becomes  the  total  energy.  The  error  introduced  by  the  numerical  integrator  tn  the  integral  of  motion  during 
consecutive  time  steps  provides  a  good  measure  of  the  local  integration  error,  and  permits  a  simple  and  reliable  strategy 
for  varying  the  time  step.  It  is  also  shown  how  the  energy  stored  in  the  penalty  system  is  suitable  to  measure  the  local 
integration  error  Overall,  the  method  is  efficient  and  powerful;  it  is  suitable  for  stiff  and  non-stiff  systems,  robust  for 
all  time  step  sizes,  and  it  works  for  singular  configurations,  redundant  constraints  and  topology  changes.  Also  the 
constraints  in  positions,  velocities  and  accelerations  are  satisfied  during  the  simulation  process.  The  method  is  robust 
in  the  sense  that  it  becomes  more  accurate  as  the  time  step  size  decreases. 

A  section  is  devoted  at  the  end  of  the  paper  to  present  numerical  simulations  that  illustrate  the  performance  of  the 
proposed  method. 


1  INTRODUCTION 

Kinematics  and  dynamic  simulations  of  multibody  systems  allow  the  accurate  prediction  of  the  behavior  of  heavy 
machinery,  spacecraft,  automobile  suspensions  and  steering  systems,  graphic  arts  and  textile  machinery,  robots, 
packaging  machinery,  machine  tools,  etc.  The  first  issue  to  consider  in  the  simulation  process  is  that  of 
modeling  the  system;  that  is,  the  selection  of  a  set  of  parameters  or  coordinates  that  will  allow  to  unequivocally 
define  at  all  times  the  position,  velocity  and  acceleration  of  the  system.  The  most  useful  kinds  of  coordinates 
currently  used  to  define  the  motion  of  multibody  systems  are  relative  coordinates,  reference  point  (or  Cartesian) 
coordinates,  and  natural  (or  fully  Cartesian)  coordinates.  These  coordinates,  when  combined  with  the  principles  of 
dynamics,  lead  to  the  final  form  of  the  equations  of  motion.  Dynamic  principles  such  as  Lagrange’s  formulation, 
Newton’s  Laws,  canonical  equations  of  Hamilton,  Virtual  Power,  Hamilton’s  Principle  and  Gibbs-Apell 
equations,  constitute  the  basis  for  the  formulations  of  multibody  dynamics1-2'3'4.  The  choice  of  dynamic 
formulation  determines  the  subsequent  choice  of  numerical  integration  schemes. 

The  method  of  Lagrange’s  multipliers  leads  to  a  representation  of  the  equations  of  motion  in  descriptor  form 
constituting  a  set  of  index-3*  differential  algebraic  equations  (DAE).  The  addition  of  stabilization  techniques,  such 
as  the  method  of  Baumgarte5,  reduces  the  index  and  makes  the  solution  tractable  by  means  of  standard  ordinary 
differential  equations  (ODE)  solvers,  however,  it  does  not  provide  full  constraint  satisfaction,  leads  to  a  limited 
control  of  accuracy,  and  in  addition  provides  no  way  for  choosing  the  values  of  the  coefficients  used  by  the 
method.  An  augmented  Lagrangian  formulation  with  projections6  has  been  proposed  which,  in  addition  to 
transforming  the  set  of  equations  into  a  stabilized  set,  is  solvable  by  standard  ODE  methods  and  assure  Lyapunov 
stability  of  the  simulation  process7.  This  method  also  has  the  advantages  of  being  robust  under  singular 
configurations,  topology  changes  and  with  redundant  constraints,  and  provides  full  constraint  satisfaction. 


'  Following  the  notation  used  by  Brenan  et  al.  (1989).  given  a  DAE  F(f.y.y')=0.  we  call  mdex  of  that  DAE  the  minimum  number  of 
times  that  all  or  part  of  the  original  DAE  must  be  differentiated  with  respect  to  the  independent  variable  (in  this  case  ()  in  order  to 
determine  the  derivative  of  the  function.  y\  as  a  continuous  function  of  y  and  t 


1 


State-space  methods,  such  as  coordinate  partitioning",  Kane  s  method  and  virtual  power  with  piojection 
matrices10,  transform  the  equations  of  motion  to  a  minimum  set  of  coordinates  that  are  Greedy  by 

methods.  State-space  methods  may  also  be  obtained  by  means  of  velocity  transformations  ■  State-space 

methods' are  more  suitable  for  ODE  integration  than  the  descriptor  counterparts  at  the  expense  of  solving  the 
velocity  and  position  problem  at  each  time  step.  However,  they  do  not  handle  well  either  topology  changes  or 
singular  configurations.  In  addition,  obtaining  the  tangent  matrix  required  for  .stiff  integrators  is  extremely 
involved  and  practically  unfeasible  in  analytical  form  using  state-space  methods.  Consequently,  they  will  not  be 
a  good  choice  when  the  system  has  numerical  stiffness.  Typically,  numerical  stiffness  in  multibody  dynamics 
anses  from  the  modeling  of  leaf  springs,  shock  absorbers,  bushings,  contact  and  impact  problems. 

The  numerical  mathematics  community  has  sought  solutions  to  the  index  reduction  problem  and  has 
proposed  many  different  ways.  Recent  advances  have  been  made  which  have  yielded  stable  index  reduction 
methods  and  accurate  ways  of  projecting  the  DAE  onto  the  underlying  ODE  for  more  stable  and  accurate 
solutions.  Key  developments  are  the  work  of  Brenan,  et  ah'6,  Griepentrog,  et  ah'  ,  Fuhrer  and  Leimkuhler  , 
Hairer  and  Wanner19  (developers  of  RADAU5),  Lubick20'2'  (developer  of  MEXX),  Petzold22  (developer  of  DASSL) 
and  Arnold23'24  (developer  of  HEX5). 

In  regard  to  numerical  integration,  the  backward  difference  formula  (BDF)  methods  have  been  customarily 
used  for  the  solution  of  differential  algebraic  equations  because  the  artificial  damping  thereby  introduced  helps  to 
stabilize  the  solution  and  provides  convergence  particularly  in  the  index-3  setting.  However,  the  actual 
implementation  of  BDF  algorithms  in  general-purpose  solvers  is  not  free  from  numerical  difficulties,  which 
become  more  acute  for  index-3  when  the  time  step  size  is  smaller  than  10'4-10  5  seconds.  The  main  difficulties 
are: 

•  For  an  index-m  DAE  the  tangent  or  quasi-tangent  matrix  used  in  the  Newton-Raphson  iteration  has  a 
condition  number  of  order  0{MAr)'<'.  Hence,  the  practical  implementation  of  the  method  is  bound  to  have 


large  round-off  errors  for  small  lime  steps  (usually  starting  at  At~\Os). 

•  Instabilities  may  result  from  sudden  changes  in  system  variables  and  constraints,  such  as  impacts,  sudden 
appearances  or  disappearances  of  constraints  and  topology  changes.  Any  time  there  is  a  discontinuity  in  the 
response  the  multi-step  BDF  tries  to  fit  a  polynomial  through  the  discontinuity  and  therefore  the  time  step 
size  must  be  severely  reduced.  As  explained  in  the  previous  point,  this  results  in  an  ill-conditioned  iteration 
matrix.  Consequently  the  Newton-Raphson  iteration  may  end  up  near  a  solution  and  yet  not  be  able  to 
converge  to  it.  These  problems  can  be  circumvented,  but  at  the  expense  of  re-initializing  the  integration,  thus 
producing  delays  in  the  integration  process. 

Due  to  the  reasons  suggested  above  we  propose  a  multi-index,  augmented  Lagrangian  formulation  of  the 
equations  of  motion  in  descriptor  form  for  index- 1  and  index-3.  The  index-3  form  is  more  efficient  than  the  index- 
1;  however,  for  time  step  sizes  smaller  than  10'5  the  ill-conditioning  of  the  tangent  matrix  in  index-3  affects  the 
performance  of  the  method,  while  the  index- 1  form  becomes  more  accurate  and  robust25.  A  single-step  numerical 
integration  scheme  with  variable  index  and  time  step  size  is  developed  based  on  the  strategy  mentioned  above. 
Mass-orthogonal  projections  to  the  constraint  space,  as  described  in  the  work  of  Bayo  and  Ledesma  ,  are 
performed  to  assure  constraint  satisfaction  to  machine  precision  during  the  integration  process.  At  the  end  of  the 
paper,  numerical  simulations  are  presented  that  illustrate  the  performance  of  the  proposed  approach. 


2.  PRELIMINARIES  ON  MULTI-BODY  DYNAMIC  ANALYSIS  IN  DESCRIPTOR  FORM 


2. 1  Formulation  in  Fully  Cartesian  coordinates 

Let  us  consider  a  multibody  system  whose  configuration  is  characterized  by  n  fully  Cartesian  (or  natural) 
coordinates2  denoted  by  vector  q  that  are  interrelated  through  the  m  holonomic  kinematic  constraint  conditions: 

O(q.r)  =  0 

The  use  of  the  Principle  of  Virtual  Power  directly  leads  to  the  equations  of  motion: 

<5qr(Mq-Q(q,q)  +  <t>qrX.)  =  0  (2> 

which  for  a  general  multibody  system  leads  to: 

Mq  +  OqrX  =  Q(q.q)  (3) 

2 


where  M  is  the  mass  matrix;  Q  is  the  vector  of  the  external,  conservative  and  non-conservative  forces;  Oq  ,s  the 
Jacobian  of  the  constraint  equations,  and  X  is  a  vector  containing  the  Lagrange's  multipliers.  The  use  of  fully 
Cartesian  coordinates  leads  to  a  constant  mass  matnx  M,  and  the  absence. of. veloc.ty  dependent  inert*  forces  m 

die  vector  Q.  r 

Equations  (1)  and  (3)  constitute  a  set  of  n  +  m  mixed  DAE's  of  index  three16,  with  q  and  X  as  unknowns.  It 
is  a  common  practice  in  multibody  dynamics  to  differentiate  twice  the  constraints,. thus  transforming  the  equation 
to  index- 1,  and  append  the  resulting  equations  to  (3)  to  yield; 


M  < 
<X>  0 


-4>q-Of 


(4) 


These  equations  can  now  be  integrated  using  standard  numerical  integration  techniques  with  each  function 
evaluation  performed  using  equation  (4).  In  addition,  equation  (4)  may  also  be  easily  modified  to  include 
Baumgarte  stabilization1. 


2.2  Index-1  Augmented  Lagrangian  Formulation  with  Prelections 

Recently,  a  method  has  been  presented  by  Bayo  and  Ledesma6  which  uses  an  index- 1  augmented  Lagrangian 
method  with  mass-orthogonal  projections  of  the  positions  and  velocities  to  their  constraint  manifolds.  This 
formulation  leads  to  the  following  equations  of  motion; 

Mq  +  4>qro4>  +  4>q7V  =  Q(q,q)  (5) 

where  X’  are  the  Lagrange  multipliers.  Introducing  the  expression  4>  =  <t>qq  +  Oqq  +  O,  in  equation  (5)  the 
following  equation  is  obtained; 

[M  +  4>^o4>q  ]q  +  <3>qrX*  =  Q(q,  q)  -  4^a(4>qq  +  )  (6) 

It  is  important  to  note  that  there  is  a  substantial  difference  between  equation  (6)  and  the  Lagrange  s 
multiplier  approach  represented  by  equation  (4).  The  leading  matrix  of  equation  (4)  becomes  singular  in  singular 
configurations,  changing  topologies  and  in  the  presence  of  redundant  constraints.  However,  although  the  mass 
matrix  M  is  in  general  positive  semi-definite,  the  leading  matrix  of  equation  (6),  [M  +  4>qra4)q  ],  is  always 
positive  definite,  which  means  that  it  can  always  be  factored,  even  in  singular  positions,  changing  topologies 
and/or  with  redundant  constraints26. 

The  augmented  Lagrangian  method  comprises  a  combination  of  penalty  terms  and  Lagrange  multipliers,  and 
they  require  an  iteration  procedure  for  the  solution  process27.  The  iteration  yields  the  unknown  multipliers  X  ,  as 
follows: 

x;+1=x;+a<i>+1,  1  =  0, 1,2,..  (7) 


with  3^=0  for  the  first  iteration.  Equation  (7)  physically  represents  the  introduction  at  iteration  i  +  1  of  forces 
that  tend  to  compensate  for  the  addition  of  all  the  constraint  terms  which  are  not  exactly  zero.  Experience  shows 
that  when  the  constraints  are  scaled  to  unity,  penalty  factors  ranging  from  105  to  107  give  good  convergence 
rates,  and  only  1  to  2  iterations  are  required  to  converge  to  the  solution.  The  leading  matrix  remains  constant 
during  the  iteration  process. 

As  a  result  of  using  the  index-1  formulation,  the  solution  of  equation  (6)  yields  a  set  of  accelerations  that 
not  only  satisfies  dynamic  equilibrium  but  also  the  constraint  conditions  <J>  =  0.  As  a  consequence,  a  projection 
in  q  to  satisfy  the  corresponding  acceleration  constraint  conditions  is  not  necessary.  On  the  other  hand,  the 
constraints  4>  =  0  and  <t>  =  0  may  not  be  satisfied  since  they  have  not  been  directly  introduced  in  the 
formulation.  Consequently,  a  projection  can  be  done  to  obtain  a  new  set  of  q  and  q  that  will  satisfy  <t>  =  0  and 
<i>  =  0,  respectively.  The  projections  in  q  and  q  may  be  done  as  indicated  right  below. 


2.2. 1  Projection  in  q 

In  order  to  satisfy  the  constraint  4>  =  0  at  each  time  step,  a  mass-orthogonal  projection  of  the  solution 
coordinates  q‘  (obtained  from  the  numerical  integration  scheme)  on  the  constraint  space  is  performed  In  this 
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way,  the  new  set  of  positions  q  that  are  obtained,  satisfies  the  constraint,  0-0.  This  can  be  achieved  by  the 
solution  of  lire  following  iterative  procedure6: 

(M  +  OqraOq)/lq('4,)  =-yJ')  (8) 

with 

Y (p  =  M(q(,)  -  q‘)  +  OqrX(,)  (9) 


and 


q(,+1)  =  q()  +  /dq(,+1),  X(,+l)  =  X{,)  +  oO(,+,)  (10) 

where  the  superscript  indicates  the  iteration  number.  liquations  (8)  through  (10)  can  be  used  iteratively  until 
fldqfl  <  £ ,  where  £  is  a  user-specified  tolerance. 

Note  the  important  feature  that  the  tangent  matrix  in  equation  (8)  is  identical  to  that  used  in  equation  (6)  for 
the  dynamic  analysis.  Furthermore,  since  the  solution  q  is  close  to  the  initial  values  q’,  the  projection  problem 
can  be  solved  using  a  modified  Newton-Raphson  method  with  no  need  for  updating  the  tangent  matrix  (usually 
convergence  is  achieved  in  just  one  iteration). 

2.2.2  Projection  in  q 

Similarly,  at  each  time  step  a  mass-orthogonal  projection  of  the  velocities  q*  on  the  velocity  constraint  space  is 
performed  to  obtain  a  new  set  of  velocities  q  that  satisfies  <f>=  0.  This  can  be  achieved  by  the  solution  of  the 
following  equation6: 

[m  +  ]q  =  Mq*  —  O7  (a<t>  +  a)  (11) 

where  the  Lagrange  multipliers  of  the  projection  problem  are  updated  as: 

a(i+,)  =  a0)  +  ex i>(,+1)  (12) 

A  more  efficient  implementation  of  the  projection  of  velocities  is  given  by  the  following  recursive  set  of 
equations: 

[M  +  07aOq]q(,+1)  =  Mq(,)  -  (13) 

where  the  iteration  is  started  by  setting, 

q^°*=q*  and  <^0)  =  ad>^0)  (14) 

Since  the  leading  matrix  in  equation  ( 1 3)  is  the  same  as  that  of  equation  (8)  used  for  the  projection  in  q, 
the  velocity  projections  will  require  only  successive  forward-reductions  and  back-substitutions. 

2.3  Index-3  Augmented  Ixigrangian  Formulation  with  Projections 

The  index-3  augmented  Lagrangian  formulation  uses  mass-orthogonal  projections  on  the  velocities  and 
accelerations  to  their  constraint  spaces.  This  method  leads  to  the  following  equations  of  motion6: 

Mq  +  OqraO  +  «I>qrr  =Q(q,q)  (15) 

where  X.'  are  the  Lagrange  multipliers.  The  following  iteration  process  yields  the  unknown  multipliers  X’: 

x;+1  =x;+a4i+I,  I  =  0.1.2....  (16) 

with  X‘0  =0  for  the  first  iteration.  Similar  to  the  index-1  formulation,  the  value  of  the  penalty  factor  a  affects 
the  convergence  rate.  Experience  shows  that,  when  the  constraints  are  scaled  to  unity,  penalty  factors  ranging 
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r  ,  j  7  fMrrmilaiion  and  once  again  only  I  to  2  iterations  are 
from  107  to  108  give  good  convergence  rates  for  the  index-3  formulatio  ,  S 

r'qUM  fotmulation.  .he  sole, ion  of  equation  (15)  y.cids  a  «,  of  O- .  -ha.  „0.  only 

satisfies  dynamic  equilibrium  bo,  also  the  comltr.ini  condihons  Asa  p"rojectio„ 

not  necessary,  and  only  q  and  q  need  to  be  projecte  to  sa  1  y  follows 

of  q  may  be  done  as  explained  above  (equatrons  (11)  to  (14)).  The  projecnon  of  q  is  done  as  fo 

2.3.1  Projections  in  q 

Similar  ,0  the  velocity  analysis,  the  projection  of  the  generalized  accelerattons  on  the  constratn,  manifold  can  be 
obtained  through  the  solution  of  the  following  equation6: 

[ M  +  <C*X>q  ]q  =  Mq ‘  -  < {<*{% Q  +  <*>,  1 +  K}  ( 1  ?) 

where  the  associated  Lagrange  multipliers  are  updated  by  the  following  formula: 

K('+,)  =  K(')+a<l>('+l)  (18) 

Numencal  implementation  of  the  augmented  Lagrangian  formulation  for  the  projection  of  accelerations  is 
suggested  by  the  recursive  set  of  equations: 

[M  +  <aOq]q('  +  l)  =  Mq(,)  -  <a{<t>qq+ <t>, }  (19) 

where  the  recursion  is  started  by  setting: 


;(<» 


:  q’  and  tc 


-<°)  -  c^>(0) 


(20) 


Equation  (19)  constitutes  a  system  of  algebraic  equations  in  q  which  is  solved  iteratively  until  convergen 

of  the  generalized  accelerations  is  achieved.  Since  the  leading  matrix  in  aquation  (19)  is  the  same  as  ^  a( :  o 
equation  (13)  used  for  the  velocity  projection,  the  accelerauon  analysis  will  require  only 
forward-reductions  and  back-substitutions. 

3.  NUMERICAL  SOLUTION  OF  THE  EQUATIONS  OF  MOTION 

During  the  integration  process  both  formulation  schemes,  index- 1  and  index-3,  lend  to  a  nonlinear  system  of 
simultaneous  equations  for  every  time  step  (equations  (6)  and  (15))  with  q  and  X  as  pnmtny  unknowns.  The 
solution  of  this  set  of  equations  requ, res  an  iterance  process,  commonly  based  on  a  Newton-Raphson  prKedme 
as  ,1  will  be  shown  in  the  following  paragraphs.  This  procedure  is  developed  by  combining  (6)  and  (1 5)  with  th 
difference  equations  of  the  numencal  integration  scheme. 

The  difference  equations  of  the  numerical  integration  scheme  may  be  expressed  as. 


=  a9n  +  l  _CL-I  and  9-tel  ~  ^ffn+l  9 n  +1 


(21) 


where  a  and  b  are  constants  that  depend  on  the  numerical  integrator  and  the  time  step,  qn  +  ,  and  qn+I  are  known 
quantities  that  depend  on  the  positions,  velocities  and  accelerations  at  step  n  and/or  previous  steps.  Note  that  i 
the  method  is  explicit,  a  and  b  will  be  equal  to  zero. 

3. I  Index- 1  Augmented  Lagrangian  Formulation 

Substituting  equation  (21)  in  the  equation  of  motion  (6)  for  index-1,  yields: 

[M  +  ^O^q)bqn  +  l+^'  +  +  (22) 

which  constitutes  a  set  of  n  equations  of  the  form  f(q„  +  l )  =  0  to  be  solved  for  the  unknown  posit.on  qn+1- 
Note  that,  by  introducing  the  difference  equations  (21).  the  positions  at  step  n  +  1  have  become  the  primary 
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variables  as  compared  wUh  .he  paginal  index- 1  formola.ion  <6,  which  .he  acc.lerahops  were  ,he  pnman, 

“k:::si  Ne„.„„-R»phS0„  procedn.c  can  he  applied  ,o  evalua.e  solohon  of  <2,,  by  means  of  .he  i.e.ahve 

process: 


ill  V'+1,=-f(q(,')) 
3qJ  V 


where  the  function  f  is  defined  as: 


f(q„+, }  =  Mqn  +  |  +  <  (a*  +  X*)  -  Qn+. 


and  the  tangent  matrix  is  approximated  by  the  following  matrix: 

—  =  fc(M  +  Oqr0<t>q)-Qq-aQq  (25) 

9q 

where  the  last  two  terms  of  (25)  represent  the  contribution  to  the  tangent  matrix  coming  from  the  elastic  forces 
(i  e  springs)  and  those  dependent  on  the  velocity  terms  (i.e.  dash-pots).  In  order  to  obtain  the  tangent  ma  nx 
equation  (25),  higher  order  terms  like-products  of  4>qq  with  <f>  have  been  neglected.  These  terms  are  negl.gib 

compared  to  those  of  equation  (24)  (see  Garcia  de  Jalon  and  Bayo\  Chapter  8). 

3. 2  Index-3  Augmented  Lagrangian  Formulation 

The  substitution  of  equation  (21)  into  equation  (15)  for  index-3  yields 

fcMqn+1+<(c^  +  ^)  =  Qn+,+Mqn  (26) 

which  again  constitutes  a  set  of  nonlinear  algebraic  equations  with  qn+1  as  primary  unknowns.  A  modified 
Newton  Raphson  iteration  is  used  again  to  evaluate  the  solution  by  means  of  the  iterative  process  outlined  above 
with  the  function  f  defined  as: 

f(q„+l)  =  Mq„+1  +  <  (<*!>  +  X* )  -  Q„+i  (27) 

The  reader  should  notice  that  equations  (24)  and  (27)  are  not  the  same:  the  constraints  <J>  are  doubly  differentiated 
in  equation  (24)  whereas  they  are  not  in  (27).  The  tangent  matrix  is  approximated  by 

—  =  fcM  +  d>qroOq-Qq  -aQq  (28) 

3q 


3.3  Numerical  behavior  of  the  formulations.  A  Simple  Numerical  Example. 

In  order  to  illustrate  the  behavior  of  the  methods  outlined  above  for  different  time  steps,  we  perform i  the 
simulation  of  a  double  pendulum  which  moves  from  at  rest  conditions  in  the  horizontal  position  under  gr  y 
effects  (Figure  1).  Each  link  of  the  pendulum  has  a  distributed  unit  mass  and  a  unit  lengt  •  e  usc 
fmd.na.5  fo.  L  modeling  p.ocess’  and  for  fee  .m.gra.ion  we  «Se  .he  bapepoida!  ru,«  wh.ch  ,s te  « £ 
second  order,  implicit,  A-stable  method  and  it  is  also  energy  preserving  in  the  linear  regime.  The  total 
simulation  is  10  seconds  and  the  penalty  factor  is  107. 

Table  1  shows  the  maximum  error  in  the  energy  norm  resulting  from  each  method,  index- 1  and  index-3,  tor 
different  time  steps  of  integration.  It  is  important  to  point  out  that  the  index-3  formulation  does  not  converge 
time  steps  enual  or  smaller  to  lO’5  seconds  when  no  pivoting  is  used.  This  fact  is  due  as  mentioned  , 
ill-conditioning  of  the  tangent  matrix  (equation  (28)).  The  problem  may  be  partially  circumvented  by  either 
increasing  the  penally  factor  or  by  using  a  pivoting  strategy  in  the  tnangulanzation  of  the  tangent  matrix. 

Table  1  shows  how  the  pivoting  helps  the  Convergence,  and  extends  the  range  in  which  the  index-3  me 
can  be  applied,  in  this  case  up  to  10'5  seconds.  However,  for  time  step  of  I06  seconds  (despite  the  time  taking 
the  integration)  the  error  starts  increasing  even  with  pivoting,  and  the  method  lacks  robustness. 


Conversely,  the  index-1  method  is  more  robust  than  the  index-3  for  time  steps  size  smaller  nan  10 
seconds.  Robustness  is  meant  in  the  sense  that  it  increases  in  accuracy,  convergence  and  does  no  show 
numerical  ill-conditioning  as  the  time  step  decreases.  On  the  other  hand,  s.nce  the  index-3  formula  .on  is  simpler 
and  involves  a  smaller  amount  of  nonlinear  terms  is  more  efficient  numerically  than  the  index-!  (see  the 
execution  time  columns  in  Table  1 ).  This  trend  is  even  more  pronounced  for  large  scale  systems  . 

It  may  be  concluded,  therefore,  that  a  multi-index  formulation  would  take  advantage  of  the  good  qualities  of 
both  methods  for  different  time  step  sizes.  Index-3  behaves  well  in  accuracy  and  efficiency,  and  .s  the  right  choice 
for  large  lime  steps;  however,  as  the  size  of  the  time  step  is  decreased,  starting  in  the  neighborhood  of  4/- 10 
seconds,  the  index-1  method  becomes  more  accurate  and  robust.  Aswitch  mechanism  can  be  implemented  within 
the  framework  of  a  variable  time  stepping  algorithm,  which  will  be  described  next. 


Time  step  size 

Index- 1 

Index-3  Without  Pivoting 

Index-3  Wi 

th  Pivoting 

Error 

Exec,  time 

Error 

Exec,  time 

Error 

Exec,  time 

At  =  10'2 

3.8e-l 

0.4  s 

2. 1  e- 1 

0.4  s 

2.2e-l 

0.4  s 

Ar=  10'3 

4.68e-3 

3.14  s 

4.6e-3 

1 .6  s 

4.8e-3 

1.67  s 

L- 

II 

p 

4.3e-5 

26.6  s 

2.0e-5 

15.2  s 

4.4e-5 

15.9  s 

o 

II 

-5 

1.25e-7 

254.6  s 

Wrong  res. 

- 

3.6e-7 

92.5  s 

At  =  10'6 

1.0e-8 

2550 

Wrong  res. 

- 

..  1.1  e-4 

922.9  s 

Table  I.  Comparative  results  for  the  double  pendulum  using  Index- 1  and  Index-3. 


4.  A  MULTI-INDEX  VARIABLE  SINGLE  STEP  METHOD 


As  mentioned  above  the  index- 1  method  is  more  accurate  than  the  index-3  for  small  time  steps,  and  the  index-3 
method  is  more  efficient  than  the  index- 1  (with  same  accuracy)  for  large  time  steps.  Therefore,  a  multi-index 
formulation  would  lake  advantage  of  the  good  qualities  of  both  methods  for  different  time  step  sizes.  Index-3  is 
the  right  choice  for  time  steps  between  I  O'2  and  104  (where  it  is  accurate  and  efficient);  however,  as  we  decrease 
the  size  of  the  time  step,  starting  in  the  neighborhood  of  At  -  10^(with  no  pivoting)  or  At  =  10  (with 
pivoting),  the  index- 1  method  becomes  more  accurate  and  robust. 

When  performing  a  simulation,  it  is  desirable  to  obtain  the  most  accurate  results  with  the  least 
computational  effort.  Therefore,  it  is  preferable  to  use  an  index-3  approach  rather  than  an  index- 1  whenever 
possible  {possible  in  this  case  means  when  the  accuracy  in  the  results  is  sufficient  for  the  given  purposes). 
Hence,  choosing  an  integration  time  step  is  a  critical  task:  if  the  At  is  too  small,  the  integration  process  will  last 
more  than  needed,  probably  with  no  evident  benefits  in  the  accuracy  of  the  numerical  results.  Conversely,  if  the 
time  step  is  tod  large,  the  integration  process  maydiverge  or  the  results  obtained  may  be  wrong.  Consequently, 
the  optimum  combination  of  integration  time  step  and  index  approach  for  each  problem  should  be  sought. 

Considering  a  single  and  fixed  time  step  algorithm,  the  best  combination  of  index  approach  and  time  step 
will  be  determined  by  the  nature  of  the  problem  and  by  the  worst  function  conditioning  in  the  sense  of  its  time 
history.  That  is,  if  the  function  being  integrated  is  smooth  all  over  the  time  interval,  a  moderate  time  step  size 
combined  with  an  index-3  approach  will  be  a  good  choice.  However,  in  cases  where  the  function  has  sharp  zones, 
even  of  short  duration,  the  whole  integration  process  would  have  to  be  carried  out  using  a  small  time  step, 
otherwise  it  would  fail.  The  choice  for  the  right  index  would  be  determined  by  the  step  size.  As  a  general  rule  if 
the  time  step  size  is  smaller  than  10  '  seconds  an  index- 1  should  be  the  choice.  Conversely,  for  time  step  sizes 
over  10'4  seconds,  the  choice  should  be  index-3  There  is  a  gap  between  10'4  and  10  s  in  which  both  methods 
perform  well  and  the  choice  becomes  irrelevant 

A  variable  time  step  integrator  would  certainly  lead  to  an  economy  in  CPU  time  for  the  case  considered 
above.  If  the  time  step  is  variable,  then  it  will  be  possible  to  accommodate  it  to  obtain  constant  local  error  or  to 
maintain  this  error  below  some  given  reference  value.  Classical  approaches  in  ODE  establish  a  variable  time  step 
strategy  based  on  a  measure  of  the  local  truncation  error,  which  is  determined  by  evaluating  either  the  integrated 
functions  using  two  different  order  methods  (the  Runge-Kutta-Fehlberg,  for  example),  or  by  integrating  two 
successive  time  intervals  with  step  size  h„  and  hjl.  Any  of  these  strategies  require  many  more  extra  function 
evaluations,  thus  compromising  the  numerical  efficiency  of  the  method. 

A  first  thought  is  to  use  the  kinetic  energy  stored  in  the  penalty  system  as  a  measure  of  the  local  integration 
error  This  energy  is  evaluated  in  the  following  way 
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(29) 


K  =  —  a <T>  <£> 
°  2 


which  can  be  computed  easily  and  efficiently.  Figures  2,  3  and  4  show  the  vertical  acceleration  of  the  end  point 
in  the  double  bar  pendulum,  the  penalty  system  kinetic  energy,  and  the  total  system  energy,  respective  y.  t  may 
be  observed  that  the  acceleration  is  very  spiky  because  the  second  bar  of  the  double  pendulum  undergoes  a  very 
strong  motion.  Also,  when  sudden  changes  in  the  acceleration  take  place  the  penalty  kinetic  energy  increases  and 
the  total  system  energy  varies.  Hence,  qualitatively  the  energy  in  the  penalty  system  can  be  taken  as  a  measure  of 
the  local  integration  error,  provided  it  increases  when  the  system  encounters  more  difficulties  in  satisfying  the 


constraints. 

Although  good  in  principle,  this  manner  of  evaluating  the  local  integration  error  has  an  important  drawback: 
it  is  quite  difficult  to  establish  a  general  quantitative  relationship  between  the  local  integration  error  and  the 
energy  stored  in  the  penalty  system. 

Due  to  this  reason,  and  taking  into  account  the  existing  correlation  between  the  energy  and  the  local  error 
(see  Figure  4  and  Figure  5),  we  propose  to  use  an  integral  of  motion  (system  total  energy  in  conservative 
systems)  as  a  measure  of  the  local  integration  error.  For  non-conservative  systems  and  considering  the  use  of 
fully  Cartesian  coordinates,  the  following  integral  of  motion  (energy  invariant)  may  be  established, 
premultiplying  the  equations  of  motion  by  qT  and  integrating  over  time  the  following  expression  may  be 


obtained: 


J,qr(Mq  +  d>qrX-Q)dr  =  constant  (30) 

Since  the  constraint  forces  are  composed  of  sets  of  internal  forces  that  do  not  produce  work,  the  integration  of 
this  equation  leads  to: 

n  =  T(t)+  v{t)-  J^qrQdr  =  constant  (31) 


Equation  (3 1 )  clearly  yields  T  +  V  =  constant  for  conservative  systems. 

5.  THE  MULTI-INDEX  AUTO-TIME-STEPPING  ALGORITHM 

Using  the  energy  invariant  described  in  the  previous  section  (equation  (31))  as  a  measure  of  the  local  integration 
error,  a  variable  time  step  size  strategy  is  proposed.  The  key  idea  is  to  modify  the  step  size  when  the  change  in 
the  energy  invariant  exceeds  an  allowable  value.  When_this  occurs  a  new  time  step  size  is  calculated  in  order  to 
achieve  the  desired  accuracy.  The  new  time  step  size  At„  can  be  expressed  as  a  function  of  the  previous  time  step 
size  At„  in  the  following  way28: 

Ain  =  rirAin  (32) 


where  qr,  is  a  parameter  defined  as: 


(33) 


where  T,  is  a  function  of  the  local  error  estimation;  £  is  the  allowed  local  error;  r  the  order  of  the  integrator  (it 
should  be  set  to  2  for  the  trapezoidal  rule)  and  v  a  safety  factor  typically  taken  as  0.8  In  fact,  good  results  are 
obtained  for  values  of  the  parameter  r  in  the  range  0.25  <  l/r  <  0.5  . 

The  strategy  to  perform  the  change  of  index  is  based  on  the  general  consideration  that  when  integrating  in 
index-3,  the  larger  the  time  step  size,  the  better  the  performance,  and  that  the  time  step  should  never  be  smaller 
than  I  O'4.  On  the  contrary,  for  the  index- 1  approach,  the  smaller  the  step  size,  the  better  the  results.  There  is  an 
overlapping  zone  ranging  from  At  =  10  4  to  10'5  in  which  both  approaches  are  suitable  In  this  zone,  it  is 
preferable  to  use  the  more  economical  index-3  rather  than  the  index- 1  approach. 

Figure  6  illustrates  in  a  flowchart  manner  the  general  strategy  established  for  the  modification  of  the  time 
step  size.  The  index  change  is  also  based  on  convergence  considerations.  If  any  integration  step  does  not 
converge,  or  the  local  error  is  larger  than  a  specified  allowable  value  Ur ,  a  time  step  reduction  is  performed.  If  the 
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problem  persists,  using  the  index-3,  after  several  consecutive  time  step  reductions  (denoted  as  maxite  in  the 
flowchart  of  Figure  6),  even  for  step  sizes  larger  than  10  4  then  it  is  necessary  to  change  to  index- 1.  n  the 
contrary  when  the  local  error  is  below  some  specified  value  L,  the  step  size  is  enlarged  and  a  change  of  index  is 
performed  (from  index-1  to  index-3)  provided  the  step  size  is  larger  than  10'f  Based  on  empirical  evidence  the 
estimated  values  for  U,  and  Lr  are  around  10 5  and  10  7,  respectively. 

It  is  important  to  remark  that  prior  to  the  change  to  index- 1  from  index-3,  it  rs  necessary  to  recalculate  the 
accelerations  q  and  velocities  q  to  satisfy  equation  (5)  for  the  current  positions  q. 

6.  NUMERICAL  RESULTS 

In  this  section  two  examples  are  presented.  First  of  all  the  double  pendulum  described  earlier  in  thrs  paper  and  the 
front  suspension  system  of  an  off-road  vehicle. 

6. 1  Double  Pendulum 

As  it  has  been  written  before,  the  double  pendulum  case  consists  on  the  simulation  of  the  behavior  of  a 
double  pendulum  released  from  its  horizontal  resting  position  for  a  total  time  period  of  10  seconds.  Each  bar  of 
the  pendulum  has  a  length  of  1  meter  and  a  uniform  mass  of  1  Kg. 

•Figure  7  and  Figure  8  show  the  results  obtained  in  this  simulation.  It  may  be  seen  that  index  and  time  step 
size  changes  occur  when  sudden  variations  in  the  system  energy  invariant  take  place.  In  this  simple  example  the 
overhead  introduced  by  the  variable-time-stepping  algorithm  is  comparatively  large  to  the  total  CPU  time 
consumed.  Hence  no  benefits,  in  terms  of  CPU  time,  are  obtained.  However,  considerable  accuracy  is  obtained 
when  integrating  at  variable  time  step  size,  as  seen  in  Figure  8,  when  compared  to  fixed  time  step  size. 


6.2  Off- road  vehicle  suspension  system 

The  1/4  ton  4x4  litis  vehicle29  has  been  proposed  as  a  benchmark  problem  by  the  European  automobile  industry 
to  check  multibody  dynamic  codes. 

We  perform  the  simulations  using  the  front  suspension  of  the  vehicle  (see  Figure  9).  The  model  requires  a 
total  of  23  variables,  related  through  22  constraint  equations,  since  there  is  only  one  degree  of  freedom.  The 
characteristics  of  the  suspension  are: 

•  a  leaf  spring,  modeled  as  a  linear  spring  of  stiffness  35,900  N/m  for  deformations  smaller  than  14.5 
centimeters.  When  such  deformation  is  exceeded  the  suspension  hits  a  second  spring  of  value  107  N/m,  thus 
adding  a  much  higher  stiffness  to  the  suspension; 

•  a  shock  absorber,  which  provides  an  elastic  force  in  Newtons  due  to  an  external  polymer,  given  by, 

Fs  =  -4.0092  106  +  2.8397  107x  -  6.7061  lO7*2  +5.2796  10  V  (34) 

and  a  damping  force  also  in  Newtons  given  by  the  following  formula: 

F0  =  -208.21  +922.1 5v  v<-0.2m/s 

Fd  =4972.8 135v+  16977. 86v2  -  29916. 125v3  -  197825^  -0.2  <  v  <  0.21  m/s  (35) 

Fd  =959.5819  +  817. 3635v  v>0.21m/s 


where  the  distance  x  is  in  meters; 

•  a  tire,  modeled  by  means  of  a  linear  vertical  spring  of  stiffness  460,000  N/m. 

The  simulation  is  performed  as  follows:  the  suspension  starts  moving  from  at  rest  conditions  at  a  speed  of  5 
m/s,  and  its  position  does  not  correspond  to  the  static  equilibrium  Thus,  it  freely  oscillates  until  the  static 
equilibrium  position  is  reached.  After  three  seconds,  the  suspension  goes  over  a  road  bump,  defined  by  a  cosine 
function,  and  afterwards  it  freely  oscillates  until  the  equilibrium  is  reached  again.  The  complete  analysis  lasts  for 


6  seconds. 

This  simulation  is  carried  out  under  two  different  conditions.  The  first  one  is  done  using  a  index-3  approach 
with  a  fixed  time  step  of  10  2  seconds.  The  second  simulation  is  performed  at  variable  lime  step  and  multi-index, 
starting  with  index-3  and  a  step  size  of  10 '3  seconds.  The  tolerance  values  U,  and  Lr  are  set  to  10 6  and  10  , 
respectively.  The  first  simulation  took  8.41  seconds  of  CPU  time  and  the  second  took  36.07  seconds. 
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Figure  10  shows  the  time  history  of  the  time  step  size.  It  may  be  clearly  seen  that  for  the  specified 
tolerance  it  is  necessary  to  shorten  the  time  step  size  only  when  the  wheel  is  passing  over  tie  ump.  I  e 
beginning  the  step  size  is  quickly  adjusted  to  speed  up  the  process  while  maintaining  the  local  error  be  ween  the 

allowable  limits. 

Figure  11  and  Figure  12  show  the  total  system  energy  in  both  cases.  It  is  important  to  remark  that  only 
when  translating  the  ongin  and  scaling  the  energy  represented  in  the  vertical  axis  (Figure  12)  one  can  notice  its 
change.  TTte  biggest  deviation  from  the  initial  value  (425,128  J  vs.  100  J)  is  about  0.02%  of  the  total  system 
energy  in  both  cases. 

It  is  important  to  remark  that  similar  precision  in  total  system  energy  conservation  can  be  obtained 
performing  the  simulation  with  an  index- 1  approach  and  a  fixed  time  step  size  of  5-10  seconds.  In  that  case  the 
simulation  takes  155  seconds  (over  400%  more  than  in  the  variable  index-vanable  step  case),  thus  corroborating 
the  benefits  of  using  the  proposed  multi-index  variable  time  step  procedure. 

As  a  result  of  this  example,  it  is  clear  that  a  multi-index- 1-3  approach  with  variable  time  step  provides  an 
excellent  solution,  since  the  more  efficient  index-3  method  is  used  during  a  large  part  of  the  motion,  and  the 
index- 1  only  at  that  part  where  the  time  step  decreases  below  the  critical  limit  where  the  index-3  loses  precision 

and  accuracy. 


7.  CONCLUSIONS 

The  proposed  multi-index  formulation  shows  a  behavior  that  could  be  defined  as  complementary.  On  one  hand, 
index-3  is  very  accurate  and  more  efficient  than  index- 1  for  large  time  steps,  whereas,  index-3  provides  wrong 
results  for  very  small  time  steps  due  to  numerical  ill  conditioning.  In  order  to  accommodate  all  these  features  a 
multi-index  variable  time  step  size  strategy  has  been  devised  based  on  the  following  criteria: 

1 .  Whenever  possible  the  time  step  is  increased. 

2.  The  time  step  size  is  decreased  based  on  the  measure  of  the  total  energy  invariant,  which  has  been  showed  to 
maintain  a  relationship  with  the  local  error. 

3.  The  threshold  to  change  from  index-3  to  index- 1  is  set  to  10J. 

These  are  the  major  conclusions: 

•  The  system  energy  invariant  has  proved  to  be  a  good  measure  of  the  local  integration  error.  This  feature 
makes  unnecessary  the  use  of  traditional  (but  less  efficient)  error  criteria  used  in  standard  numerical  integration 
methods. 

•  The  variable  time  step  strategy  and  the  switch  from  one  method  to  the  other  do  not  cause  problems  during  the 
integration  process. 

•  The  use  of  the  proposed  technique  allows  a  substantial  speedup  gain  in  CPU  time  for  large  scale  systems, 
thus  helping  to  achieve  real-time  behavior. 

•  The  method  is  general  and  can  also  be  applied  to  solve  the  dynamics  of  flexible  multibodies. 
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Figure  1.  Double  pendulum  initial  position. 

Figure  2.  Vertical  acceleration  of  the  double  pendulum  end  point. 

Figure  3.  Double  pendulum  kinetic  energy  of  penalty  system. 

Figure  4.  Total  system  energy  vs.  time 
Figure  5.  Local  integration  error  vs.  time 

Figure  6.  Multi-index  and  time-stepping  strategy  flowchart  block  diag-am. 
Figure  7.  Time  history  of  the  time  step  size. 

Figure  8.  Total  system  energy. 

Figure  9.  Model  of  the  litis  suspension  system. 
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Figure  10.  Time  history  of  the  time  step  size. 

Figure  11.  Total  system  energy. 

Figure  12.  Amplification  of  the  total  system  energy. 
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Figure  1.  Double  pendulum  initial  position. 


Figure  3.  Double  pendulum  kinetic  energy  of  penalty  system 
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Figure  4.  Total  system  energy  vs.  time 


Figure  5.  Local  integration  error  vs.  time 
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Figure  6.  Multi-index  and  time-stepping  strategy  flowchart  block  diagram. 
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Figure  7.  Time  history  of  the  time  step  size. 
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Figure  8.  Total  system  energy 


Figure  9.  Model  of  the  litis  suspension  system. 
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Figure  10.  Time  history  of  the  time  step  size 
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Figure  11.  Total  system  energy. 


Figure  12.  Amplification  of  the  total  system  energy. 
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Abstract.  Real-time  dynamic  simulation  of  large,  realistic  and  complex  multibody  systems  is 
essential  to  develop  modem  technologies  such  as  virtual  prototyping,  man-in-the-loop  simulators  and 
intelligent  vehicle  control  systems.  In  order  to  achieve  real-time  performance,  current  commercial 
codes  require  the  use  of  large  costly  computers,  thus  limiting  the  number  of  potential  users. 

This  paper  shows  that  real-time  can  be  achieved  on  medium-size  workstations  if,  on  one  hand,  an 
adequate  combination  of  modeling,  dynamic  formulation  and  numerical  integration  scheme  is  selected 
and,  on  the  other  hand,  advantage  is  taken  from  sparse  matrix  technology  and  parallel  computing.  A 
study  of  space-state  and  descriptor  methods  involving  the  dynamics  of  a  whole  vehicle  model  is 
earned  out,  and,  as  conclusion,  two  methods  are  proposed  as  the  best  candidates  for  real-time 

simulation. 

Key  words:  real-time-simulation,  multibody  dynamics,  solution  algonthms,  parallel  computing, 
sparse  matrix  technologies. 


1.  Introduction  and  Background 

Simulation  tools  are  becoming  more  and  more  relevant  in  mechanical  design  as  they  allow  for  a 
reduction  of  the  design-cycle,  thus  leading  to  products  at  reduced  costs  and  with  earlier  presence  in 
the  market  place.  More  specifically,  computer  aided  dynamics  of  multibody  systems  is  of  great 
interest  for  automotive,  aerospace,  robotics.,  biomechanical  and  military  industries.  In  a  number  of 
applications,  the  time  spent  by  the  computer  in  performing  the  analysis  is  not  a  crucial  matter. 
However,  there  are  applications  that  cannot  be  developed  unless  real-time  simulation  is  achieved.  This 
group  encompasses  hardware-in-the-loop  settings,  where  a  real  component  may  be  tested;  man-m- 
the-loop  devices,  used  for  training  purposes;  virtual  prototyping,  that  allows  the  designer  to  interact 
with  his  model  and  immediately  see  the  results  of  a  what-if  analysis;  and  intelligent  vehicle  control 
systems,  that  lead  to  safer  and  more  comfortable  transportauon  vehicles.  Available  codes  and 
algorithms  are  capable  of  reaching  real-dme  performance  on  convendonal  computers  when  managing 
small  and  academic  examples,  but  large,  realistic  and  complex  models  lead  to  very  large  calculadon 
efforts  which  require  powerful  computers  and  efficient  codes  to  meet  real-time. 

In  a  previous  paper  [1],  the  authors  earned  out  a  study  on  the  different  factors  involved  in  the 
simulation  process:  modeling,  dynamic  fomiulation  and  integration  procedures;  and  concluded  that 
the  adequate  combination  of  these  factors  depends  on  the  properties  of  the  system  to  be  analyzed. 
Therefore,  it  is  not  correct  to  talk  about  the  best  method  for  all  cases  but  about  the  best  method  for  a 
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particular  problem,  and  consequently  the  concept  of  intelligent  simulation  was  introduced.  In  ( l],  four 
methods  were  developed  and  tested  to  find  out  what  kind  of  problems  they  managed  better.  For  this 
purpose,  a  benchmark  library  containing  open  and  closed  kinematic  loops,  mechanisms  with  singular 
configurations  and  stiff  systems  was  also  set  up.  Results  led  to  conclude  that  two  of  the  methods 
were  notably  superior  to  the  others:  an  index-3  augmented  lagrangian  formulation  with  mass- 
orthogonal  projections  (a  descriptor  method),  and  a  fully-recursive  formulation  based  on  relative 
coordinates  (a  state-space  method).  However,  none  of  the  examples  examined  in  [1J  was  large 
enough  to  be  able  to  generalize  the  conclusions  for  complex  systems.  Realistic  models  of  vehicles, 
satellites,  robots  or  humans  usually  exceed  150  coordinates  in  size.  Furthermore,  parallel  compuUng 
and  sparse  matrix  solvers,  whose  influence  is  noticed  when  the  problem  size  increases,  are  factors  to 

be  taken  into  account  which  had  not  been  considered. 

These  are  die  reasons  that  led  to  develop  the  work  presented  in  this  paper.  The  two  methods 
mentioned  above,  and  a  third  one,  an  index-3  classical  Lagrangian  formulation,  commonly  used  in 
commercial  codes,  are  described  in  detail.  Emphasis  is  made  on  the  implementation  aspects  which  are 
crucial,  and  constitute  a  major  improvement  with  respect  to  reference  [1],  Then,  the  dynamics  of  the 
complete  model  of  a  4x4  military  vehicle  composed  of  168  coordinates  undergoing  several  demanding 
maneuvers,  are  solved  using  the  three  methods  under  comparison.  The  influence  of  sparse  matrix 
technology  and  parallel  computing  is  also  examined.  Based  on  the  results,  it  is  concluded  that  real¬ 
time  performance  can  be  reached  on  medium-size  workstations. 

2.  Dynamic  Formulations  and  Numerical  Implementation 

In  this  section,  the  three  methods  whose  performances  are  to  be  compared  will  be  briefly  described. 
Description  of  each  method  includes  the  form  of  the  dynamic  equations  and  the  integration  scheme. 

2  i  INDEX-3  AUGMENTED  LAGRANGIAN  FORMULATION  WITH  PROJECTIONS 

The  index-3  augmented  Lagrangian  formulation  with  mass-orthogonal  projections  may  be  found  in 
[2],  The  corresponding  equations  of  motion  are  given  by 

Mq  -t-OjaO  +  oJX*  =  Q  (D 

where  M  is  the  mass  matrix,  q  are  the  accelerations,  the  Jacobian  matrix  of  the  constraint 

equations,  a  the  penalty  factor,  the  constraints  vector,  X*  the  Lagrange  multipliers  and  Q  the 
vector  of  applied  and  velocity  dependent  inertia  forces.  The  Lagrange  multipliers  are  obtained  from 
the  following  iteration  process  [2], 

a,(-+1  —X[  +  ctOj+i ,  £=0, 1,2 —  (2) 

In  [2],  the  value  X*0  =  0  was  chosen  to  start  the  iteration.  However,  more  recent  experiences  have 
-demonstrated  that  better  convergence  is  attained  when  extrapolating  X0  from  the  Lagrange  multipliers 


,nntp,  that  sub-index  n  indicates  the  time  step,  and  sub-index  i 
already  calculated  in  previous  time  steps  (note  that  suo  inaex  n 

refers  to  the  iteration  step  within  a  time  step) 

(3) 

-  '^n-i 

>n+ 1 


(^o)  ~  ^n-1 

V  >n+ 1 


As  integration  scheme,  the  implicit  single-step  trapezoidal  title  has  been  adopted.  The  corresponding 
difference  equations  in  velocities  and  accelerations  are. 

2  ~  (  2  •  ^ 
q„+l=— ^+l+q«  Wllh  qn'VAtqn+q\ 


At 


qn+i  = — 2  qrt+l+q«  Wlth  qn 
At 


(4) 


.  — q  n  +  —  qn  +  qn 

Ur  ^  3 


(5) 


Dynamic  equilibrium  can  be  established  at  time  step  n+ 1  by  mtroducing  the  difference  equations 
(4)  and  (5)  into  the  equations  of  motion  (1),  thus  yielding 


-K  Mqn+1  +  <E>T  ,  {a®n+ 1  +  )  -  Qn+1  +  Mqn  -  0 

At1  n+1 


(6) 


Ai 

For  numerical  reasons,  the  scaling  of  equation  (6)  by  a  factor  of  —  seems  to  be  convement,  thus 


yielding 


Mqn+j  +  -r®ln+l{°&n+ 1  +  ^n+1 )  "  ~ Qn+1  +  ~~ Mq«  ~  0 


(7) 


4  4 

or,  symbolically  f(qn+i)-0. 

In  order  to  obtain  the  solution  of  this  non-linear  system,  the  widely  used  iterative  Newton- 
Raphson  method  may  be  applied 


^q 


Aq,+i  =  -[f(q)]f 


(8) 


where 


2 

[f(q)]  =  ~r(Mq  +  +  OJA.*  -  Q) 


(9) 


and  the  approximated  tangent  matrix  is: 


<*(q) 

5q 


=  M  +  ^C  +  ^«o<J>q  +  K) 


(10) 
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where  C  and  K  represent  the  contribution  of  damping  and  elastic  forces  of  the  system  provided  they 

exist.  ,  • 

A  closer  look  al  the  tangent  matrix  reveals  that  Ul-condi cloning  may  appear  when  the  Ume  step 

becomes  small  It  may  be  seen  in  equation  (10)  that  K  and  the  constraint  terms  are  multiplied  by  /Sr  , 

C  by  At  and  M  is  not  affected  by  the  step  size.  As  a  consequence  when  /Sr  reaches  small  values, 

large  round  off  errors  will  occur.  In  fact,  it  has  been  demonstrated  in  [3]  that  for  an  wdex-3 

differential  algebraic  equation  the  tangent  matrix  has  a  condition  number  of  order  1/dr 

Consequently,  the  method  is  bound  to  have  round-off  errors  for  step  sizes  smaller  than  10  . 

The  procedure  explained  above  yields  a  set  of  positions  q„+i  that  not  only  satisfies  the  equations 
of  motion  (6),  but  also  the  constraint  conditions  <t>  =  0.  However,  it  is  not  expected  that  die 
corresponding  sets  of  velocities  and  accelerations  satis:  <i>  =  0  and  *  =  0 ,  because  these  conditions 
have  not  been  imposed  in  the  solution  process.  To  overcome  this  difficulty,  mass-orthogonal 
projections  in  velocities  and  accelerations  have  been  proposed  in  [2],  However,  we  propose  a 
modification  to  the  method  of  [2]  in  order  to  get  as  leading  matrix  the  same  tangent  matrix  appearing 
in  equation  (10).  As  a  consequence,  triangulanzation  is  avoided  and  projections  in  velocities  and 
accelerations  are  performed  with  just  a  forward  reduction  and  a  back  substitution. 

If  q*  and  q*  are  the  velocities  and  accelerations  obtained  after  convergence  has  been  achieved  in 
the  Newton-Raphson  iteration,  their  cleaned  counterparts  q  and  q  are  calculated  from 


_ ,  At  At 

M  +  —  C  +  -— 
2  4 


(oJad>q  +  K 


q  = 


w  At  ^  AT  v 

M  +  — C  + - K 

2  4 


At  ^ 
4  q 


(ID 


for  the  velocities,  and 


At  ^  At 

M  H - C  H - 

2  4 


(oJafcq+K 


..  At  _  At 

M  -t - C  +  — 

2  4 


K 


At4 


■*J«(*qq  +  *t)  (12) 


for  the  accelerations. 

Sparse  matrix  technology  has  been  implemented  to  solve  all  the  linear  sets  of  equations  that  arise 
when  applying  this  method.  The  fact  that  the  leading  matrix  is  symmetric  and  positive-definite  has 
been  taken  into  account  as  well. 

2.2.  INDEX-3  CLASSICAL  LAGRANGLAN  FORMULATION  WITH  PROJECTIONS 
An  equivalent  statement  may  be  developed  for  the  classical  Lagrangian  formulation  [4],  In  this  case, 
the  equations  of  motion  are 

Mq  +  Oq^.  =  Q  (I3) 


<h  =  0 


(14) 


5 


where  the  former  relation  express  the  dynamic  equilibrium  and  the  latter  one  contains  die  constraints 
among  variables. 

Again  the  trapezoidal  role  is  chosen  to  solve  the  initial  value  problem.  The  combination  of  the 
difference  equations  (4)  and  (5)  with  equations  (13)  and  (14)  leads  to:  . 


AT 


Mqn+1  +<*>!  .An+1  "Qn+1  +  0 


4n+l 
<*>«+!=  0 


(15) 

(16) 


Equations  (15)  and  (16)  are  set  at  time  step  n+ 1.  As  for  the  augmented  Lagrangian  formulation 


described  before,  these  equations  are  scaled  by  a  factor  of  j  ,  yielding 


AT 

4 


(17) 


AT 


®n+ 1=0 


(18) 


Equations  (17)  and  (18)  may  be  written  symbolically  as  f(q„+i)  -  0  and  then  solved  by  means  of 
a  Newton-Raphson  procedure  with 

A 


[«q)]  = 


AT 


Mq  +  <&J?i-Q' 


AT 


(D 


(19) 


and  the  tangent  matrix  approximated  by 


<fr(q) 

L  <*i  J 


M  +  — C  +  ^-K 

2  4 

At 2 


(20) 


Since  only  constraint  conditions  in  positions  have  been  imposed,  it  is  not  expected  that  first  and 
second  derivatives  of  the  constraints  are  satisfied  by  velocities  and  accelerations.  Similar  to  the 
augmented  Lagrangian  formulation  seen  above,  mass-orthogonal  projections  may  be  performed  in 
such  a  way  that  the  corresponding  leading  matrices  are  identical  to  the  tangent  matrix  (20).  If  again 
q*  and  q*  stand  for  the  values  to  be  cleaned,  the  projection  for  velocities  is 
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and  for  accelerations 


At  _  4r  „ 
M  +  —  C  +  — K 
2  4 
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-T^q 
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M  H - C  +  — —  K 
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(21) 


(22) 


Again,  sparse  matrix  technology  for  symmetnc  and  positive-definite  matnces  is  used  at  the  time  of 
obtaining  simulation  results. 


2.3.  TOLLY-RECURSIVE  FORMULATION 

A  topological,  fully-recursive  algorithm  of  order  0(N)  was  developed  by  Jimenez  [5]  based  on  the 
articulated  inertia  method  originally  proposed  by  Featherstone  [6].  To  represent  the  state  of  the 
mechanism,  a  set  of  relatives  coordinates  z  are  used,  as  shown  in  Figure  1.  If  the  kinematic  chain  is 
open,  these  coordinates  constitute  a  minimum  set,  equal  in  number  to  the  degrees  of  freedom,  and 
therefore  independent;  else,  they  are  dependent  and  so  related  through  a  certain  number  of  constraint 
equations.  On  the  other  hand,  each  body  of  the  mechanism  is  defined  by  the  Cartesian  velocities 

ZT={f0T  C0T}  (23) 

where  r0  stands  for  the  velocity  of  the  material  point  of  the  body  that  is  currently  situated  at  the  origin 
.  of  the  inertial  reference  frame,  and  to  is  the  angular  velocity  of  the  body. 


Figure  1.  Configuration  of  a  multibody  system  of  6  bodies  and  6  relative  coordinates. 
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Based  on  this  kinematic  assumption,  the  principles  of  dynamics  lead  to  the  following  mass  matrix 
and  vector  of  applied  forces  for  a  certain  body 


M  = 


ml 


-mrG 


[mrj  JG  -  «rGfG  J 


Q  = 


f  -  co  x  (co  x  mrG ) 

nG  -  co  x  JGG)  +  rG  x  (f  -  CO  x  (to  x  mrG )) j 


(24) 


(25) 


where  m  is  the  mass  of  the  body,  rG  is  the  3x1  vector  containing  the  position  of  die  center  of  mass  of 
the  body,  rG  is  the  3x3  dual  antisymmetric  matrix  of  rG,  JG  is  the  3x3  inertia  tensor  referred  to  the 
center  of  mass,  and,  finally,  f  and  nG  are  the  3x1  vectors  of  applied  forces  and  moments  acting  at  the 

center  of  gravity  of  the  body,  respectively. 

As  for  any  recursive  formulation,  the  accelerations  are  calculated,  without  the  need  of  solving  a  set 
of  simultaneous  equations,  by  going  through  the  kinematic  chain  in  a  forward  and  backward  manner. 
The  first  step  consists  of  a  kinematic  transmission  starting  from  the  base  towards  die  end  of  the  chain. 
The  corresponding  recursive  expressions  in  velociues  and  accelerations  are  given  by 

Zz-  =  Z(_i  +  b(z(  (26) 

Zj  —  Z  j_j  +  b/z,  +  dj  (27) 


where  b;  is  a  6x1  vector  containing  the  value  of  the  Cartesian  velociues  of  body  i,  Z;,  when  all  the 
velociues  z  are  made  zero  except  for  zf-  =  1;  and  dz  is  also  a  6x1  vector  containing  the  difference  of 
the  Cartesian  accelerauons  Z,  -  Z,_i  when  all  the  accelerauons  z  are  made  zero.  This  first  step 
remains  the  same  whether  the  kinematic  chain  is  open  or  closed.  Obviously,  the  vectors  b(  and  d, 
depend  on  the  kind  of  kinematic  pair  that  joins  body  i-1  with  body  /'.  Equations  (26)  and  (27)  form  a 
recursive  relation  between  the  velocities  and  accelerations  of  two  consecutive  bodies,  in  terms  of  the 
relative  joint  velocity  and  acceleration. 

The  second  step  is  a  condensation  of  the  inertial  and  applied  forces  starting  from  the  end  of  the 
chain  and  progressing  towards  the  base.  If  the  kinematic  chain  is  open,  the  principle  of  virtual  power 
for  an  AMink  multibody  system  may  be  formulated  as  follows 

SZ*T(m,-Z,.-Q,-)  =  0  (28) 


where  Z*  represents  the  virtual  velocities  of  body  i.  Since  they  are  not  independent,  they  cannot  be 
eliminated  in  equation  (28).  Introducing  the  recursive  relations  (26)  and  (27)  in  equation  (28)  for  the 
body  N,  one  may  obtain: 


(29) 


Vz;T(MiZl-Qi)  + 

1=1 


Z*A,T_1[(Mw.,  +  MW)ZW_i-(QN-l+Q«'MA'dN)+M«bwiwl  + 


)]  =  o 


z*JbTN[MNbNiN  +  MNZN^-{QN-MNdNt 

The  virtual  reiat.ve  velocity  z*w  appears  in  the  thtrd  to  of  equauon  (29).  Stnee  U  is  utdependent  of 
dre  remaining  virtual  velocities,  the  bracket  multiplying  it  must  be  zero.  Thts  cond.uon  allows  us  to 
obtain  the  independent  accelerations. 

_  1  r  t  v  _  _ _  \  «  T 


zN={  bTNMNbN)  [bJlQAf-MAfd^j-b^M 


N^N-l 


(30) 


Reintroducing  Zyv  in  equation  (29),  and  defining 


Kyv  =  16  “  MA,bN(b^M^b/V) 

(31) 

> 

t 

II 

2 

> 

1 

+ 

2 

(32) 

QjV_l  =  Qn-\  +  Kjv(Q/^ 

(33) 

makes  it  possible  to  rewrite  equation  (29)  in  an  analogous  form  to  equation  (28)  but  where  all  the 
terms  affecting  body  N  have  been  eliminated 


zVfMA-  -  Q,)  +  -  Q/v-l)  -  0 


N-2 

I 

i=l 


(34) 


There  are  only  N- 1  terms  in  equation  (34)  and  MN_{  is  the  articulated  inertia  of  bodies  N- 1  and  N. 

Repeating  the  same  reasoning  in  a  recursive  manner  for  any  relative  acceleration  one  obtains  the 
accelerations  at  any  joint  i: 

v-lr 


±i  =  (bj M/b,)  [bf  (q,  -  -  b J M,Z; 


-1 


(35) 


and  the  recursive  condensation  of  applied  and  inertia  forces  which  constitute  the  second  step  for  open 
chains  is  performed  by  means  of  the  following  expressions 


K|-  =  16-^7^)"^? 

Mi-1  =  M(-l  +  K,M; 

Q,--!  =  Q.--1  +  Ki(Q«  -  *id*) 


(36) 

(37) 

(38) 
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me  third  and  las,  step  for  open  chains  is  Ore  odculauon  of  the  relative  accelerates  . t  from  dre 
base  to  the  end  of  the  chain.  Tins  calculation  is  earned  out  by  recursively  ustng  equation  (35). 

Looking  at  dre  unple  recurs, ve  procedure,  it  may  be  seen  drat  dre  number  of  —  e  operate 
grows  proportionally  with  the  nunrber  of  degrees  of  freedom  of  dre  open  chant,  an  (  me  o 

The  cons, deratron  of  branches  tn  dre  kmematre  chatn  ,s  a  stnrp.e  task.  In  the  junction  bodies,  dre 
forward  recursive  computadons  (steps  1  and  3)  must  be  split  into  two  separate  procedures  that  move 
independendy  along  each  branch.  In  the  backward  recursive  computadons  (step  2)  the  two  separate 
procedures  along  each  branch  meet  at  the  junction  body  and  yield  a  single  procedure. 

As  sard  before,  steps  2  and  3  are  not  valid  for  closed  kinemadc  chains.  However,  these  two  steps 
may  be  modified  in  order  to  tackle  closed-loop  multibody  systems.  To  obtain  .1,  they  must  be 
transformed  into  open-loop  systems  through  the  cut-jowl  method,  which  eliminates  or  cuts  a  joint  m 

each  closed-loop.  This  method  is  described  in  detail  in  [5]  and  [7], 

In  this  work  the  method  shown  in  [5]  is  followed.  When  a  joint  is  removed,  the  corresponding 
constraint  forces  (d>Jx)  and  (-O^)  must  be' ’introduced  on  both  links  connected  at  the  joint,  and 

condensed  along  with  the  other  forces.  <DZ  is  the  Jacobian  matnx  of  the  constraints  with  respect  to 
the  body  variables  Z.  and  X  are  the  Lagrange  multipliers.  A  penalty  formulation  is  used,  so  that 


X  =  a(<t>  +  2^X36  +  J22d>) 


(39) 


where  a  is  the  penalty  factor.  Equation  (39)  may  be  written  as  a  function  of  the  variables  that  define 
the  two  bodies  (say  r  and  s)  connected  at  the  cut-joint, 


X  =  a(<bzZr  -<f>zZ?  +  y) 


(40) 


y  =  4>zZr  -  4>zZ,  +  2nQ(<S>zZr  -  <S>zZs)  + Q2<t>  (41) 

Now,  the  second  and  third  steps  explained  above  for  open  chains  must  be  reformulated.  If  the  joint 
between  bodies  r  and  5  has  been  removed  to  open  a  closed-loop,  application  of  the  principle  of  virtual 
power  gives  the  following  equation 

X  Z*T(m,-Z,— Q,)  + 

'='  (42) 

i*r,s 

Z;T(MrZr  -  Qr  +  <t>zX)  +  z;t(m,Z,  -  Q,  -  OzZ)  =  0 

instead  of  equation  (28)  which  is  valid  for  open-loop  systems.  Starting  from  equation  (42),  a  similar 
process  to  that  developed  before  for  open  chains  should  be  carried  out.  Expressions  will  be  obtained 
either  for  mass  and  force  condensation  (second  step),  as  for  relative  accelerations  (third  step). 
Nevertheless,  this  process  requires  a  considerable  amount  of  effort  as  it  is  very  much  problem 
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dependent,  leading  to  involved  expressions  and  tedious  calculations.  Consequently,  general 
expressions  like  (35)-(38)  cannot  be  provrded,  and  they  are  to  be  developed  for  each  partrcu  ar 

Unlike  the  prev.ously  presented  Lagrangian  fomrulatrons  wh.clr  featured  a  Newton-Raphson 
scheme,  this  fully-recursive  method  only  allows  for  a  fixed  point  iterauon  to  be  used  for  tntegrauon 
purposes.  The  reason  is  that  introducing  the  integrator  equations  and  obtaining  the  tangent  matrix 
becomes  practically  unfeasible.  If  again  the  trapezoidal  rule  is  used,  then  positions  and  velocities 
should  be  expressed  in  terms  of  accelerations. 


At 


zn+ 1  ~~7 ~Zn+\  +zr 


with 


:  Zn  +  Atzn  + 


At1 


(43) 


At .. 


zn+\  -  T^n+1  +Zfl 


with 


At .. 


■  zn  +  — z, 
n  2 


(44) 


3.  The  Example:  Modeling  Aspects 

To  compare  the  three  formulations  presented  above,  the  dynamics  of  a  complex  and  realistic  3D 
multibody  system  have  been  analyzed.  The  chosen  example  is  the  military  4x4  Bombardier  Has 
vehicle  [8],  which  was  proposed  by  the  European  automobile  industry  as  a  benchmark  problem  to 
check  multibody  dynamic  codes.  The  vehicle  features  four  identical  suspensions  (see  Figure  2)  whose 

characteristics  are: 

-  A  shock  absorber ,  which  provides  a  damping  force  given  by  the  following  nonlinear  expression: 

=  9945.627v  +  33955.72v2  -59832.25v3  -39565 l.Ov4  [N]  for -0.2<v<0.21  m/s 
Fd  =  -416.42+  1844.3v  [N]  for  v<-0.2  m/s 
Fd  =  1919.1638+  1634.727v  [N]  forv<0.21  m/s 

and  a  elastic  force  due  to  an  external  polymer,  given  by 

Fs  =  -4.0092 *  106  +  2.8397  *  107x  - 6.7061  *  107x2  +  5.2796  *  107  x3  [N] 

where  the  distance  x  is  defined  in  meters. 

-  A  leaf  spring,  modeled  as  a  linear  spring  of  35,900  N/m.  Moreover,  a  rubber  bumper  contacts  the 
wheel  at  the  leaf  spring  connection  point  after  a  vertical  motion,  from  the  equilibrium  position  of  70 
mm.  The  bumper  stiffness  is  assumed  to  be  a  constant  value  of  107  N/m,  which  is  engaged  only  after 
the  clearance  is  used  "up.  This  last  element  is  considered  in  some  simulations,  thus  adding  a  muon 
higher  stiffness  to  the  problem,  but  not  in  others.  Its  presence  will  be  conveniently  indicated  in  the 

results  below. 
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-  A  tire,  whose  radial  stiffness  is  460,000  N/m.  The  side  force  and  aligning  torque  are  considered  by 
application  of  the  tire  Calspan  model  [9]. 

Also,  there  is  an  steering  system,  which  is  kinematically  guided  during  the  maneuvers.  The  total 

mass  of  the  vehicle  is  about  1500  Kg. 

3.1.  MODELING  FOR  THE  LAGRANGIAN  FORMULATIONS 

In  order  to  perform  the  analysis  with  Lagrangian  formulations,  the  litis  vehicle  has  been  modeled 
using  fully  Cartesian  dependent  coordinates,  also  called  natural  coordmates  [4],  For  three- 
dimensional  multibody  systems,  these  coordinates  describe  the  position  of  the  whole  mechanism  by 
means  of  the  Cartesian  coordinates  of  basic  points  and  the  Cartesian  components  of  unit  vectors ,  all 
distributed  throughout  the  elements.  Each  body  of  the  system  should  have  a  sufficient  number  of 
points  and  vectors  linked  to  it,  so  that  its  motion  is  completely  defined.  Figure  2  illustrates  how  the 
chassis,  suspension  and  steering  system  of  the  vehicle  have  been  modeled. 


Figure  2.  Modeling  of  the  litis  vehicle  in  fully  Cartesian  coordinates. 

As  may  be  seen  in  Figure  2,  basic  points,  unit  vectors  and  relative  coordinates  of  distance  have 
been  used  adding  to  a  total  number  of  168.  The  system  has  11  degrees  of  freedom:  six  free-body 
motion  degrees  of  freedom  for  the  chassis,  one  for  each  suspension  and  one  more  for  the  steering 
system.  In  consequence,  the  variables  of  the  problem  are  related  through  168-11=157  constraint 
equations.  In  practice,  only  10  from  the  1 1  degrees  of  freedom  are  true  dynamic  degrees  of  freedom 
because,  as  said  before,  the  steering  motion  is  kinematically  guided.  This  guidance  is  introduced 
through  a  new  constraint  equation,  so  that  the  final  number  of  constraints  is  158. 
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3  2  MODELING  FOR  THE  RECURSIVE  FORMULATION 

Recurs, ve  formulates  requue  the  use  of  re.ative  coordmates.  They  define  dte  position  of  each  body 
i„  relation  to  the  previous  one  in  the  krnemat.c  chain,  by  using  the  relative  degrees  of  freedom  allowed 
by  the  joint  linking  those  elements.  Figure  3  illustrates  the  way  tire  chassis,  suspension  an  steering 
system  of  the  litis  vehicle  have  been  modeled  using  this  kind  of  coordinates. 


BASE 

7/// / /V 


Figure  3.  Modeling  of  the  litis  vehicle  in  relative  coordinates. 

In  this  case,  the  number  of  variables  is  43,  distributed  as  follows:  6  for  the  chassis  (three 
translations  and  three  rotations),  9  for  each  suspension  (all  of  them  rotations  at  the  different  pairs  as 
shown  in  Figure  3)  and  1  distance  for  the  steering  system.  Since  this  last  variable  is  kinematically 
guided,  the  total  number  of  variables  amounts  to  42.  Figure  3  shows  where  the  closed-loops  have 
been  cut.  Two  cuts  are  performed  for  each  suspension.  One  of  them  needs  3+3=6  constraints  to 
establish  identity  between  points  and  directions  at  the  rotational  joint  removed  (actually,  only  two 
constraints  are  needed  for  the  direction  but  using  three  is  easier),  while  the  other  only  introduces  3 
more  constraints  as  this  time  the  joint  removed  is  spherical.  Therefore,  the  total  number  of  constraints 

is  36,  with  just  32  independent. 


4.  Simulations  Results 

In  order  to  test  the  numerical  efficiency  of  the  previously  presented  formulations,  the  litis  vehicle  is 
driven  through  three  different  simulations,  which  are  described  below. 


4.1.  SIMULATIONS  DESCRIPTION 

The  first  simulation,  that  will  be  called  Simull,  consists  of  10  seconds  of  motion  with  the  vehicle 
going  up  an  inclined  ramn  and  the  down  a  series  of  stairs  at  a  speed  of  5  m/s.  The  simulation  leads  to 
the  strong  motion  with  accelerations  of  up  5g  shown  in  Figure  4. 


Figure  4.  The  Dtis  vehicle  performing  Simull  and  vertical  acceleration  of  its  center  of  mass. 

During  the  second  simulation,  Simul2,  the  vehicle  goes  off  a  ramp  and  hits  the  ground,  leadmg  to 
the  acceleration  time  history  of  Figure  5,  that  also  lasts  for  10  seconds  at  a  speed  of  5  m/s. 


Figure  5.  The  Dtis  vehicle  performing  Simul2  and  vertical  acceleration  of  its  center  of  mass. 

In  the  third  simulation,  SimuB,  the  litis  vehicle  is  driven  through  a  series  of  columns  in  the  slalom 
maneuver.  The  simulation  lasts  for  10  seconds  at  a  vehicle  speed  of  5  m/s,  leadmg  to  the  acceleration 
time  history  of  Figure  6. 
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Figure  6.  The  litis  vehicle  performing  Simul3  and  transversal  acceleration  of  its  center  of  mass 


4.2.  RESULTS  ON  ONE- PROCESSOR  COMPUTERS 

The  three  simulations  described  above  were  earned  out  on  a  SGI  Indigo2  IMPACT  with  one 
processor  R4400SC  @  200  MHz  and  2  Mb  of  secondary  cache  memory,  using  the  three  formulations 
being  compared.  Table  I  features  two  columns  for  each  formulation:  the  first  one  gives  the  smallest 
CPU  time  achieved  with  the  corresponding  method,  while  the  second  one  shows  the  time  step  size 
which  led  to  that  best  CPU  time. 

Table  I  illustrates  the  fact  that  the  augmented  Lagrangian  formulation  (with  a  penalty  factor  of  10  ) 
is  clearly  superior  to  the  classical  one  because  it  performs  notably  faster  (between  three  and  four 
times)  and  achieves  convergence  for  greater  time  step  sizes.  On  the  other  hand,  the  fully-recursive 
formulation  shows  to  be  more  efficient  than  the  augmented  Lagrangian  method  although  it  needs  to 
take  a  smaller  time  step  size  to  converge  (due  to  the  fact  that  it  uses  fixed  point  iteration).  Thus,  in 
simulations  such  as  Simul2  the  augmented  Lagrangian  method  can  work  with  larger  time  step  sizes, 
and  performs  at  the  same  level  as  the  fully-recursive  formulation.  On  the  other  hand,  when  the 
maneuver  becomes  sharper  (such  as  in  Simull)  the  augmented  Lagrangian  method  must  keep  on 
small  step  sizes  and  the  fully-recursive  formulation  takes  a  fair  advantage. 


1 

Augmented 

Lagrangian 

Classical  Lagrangian 

Fully-Recursive 

CPU  (s) 

At  (s) 

CPU  (s) 

At  (s) 

CPU  (s) 

At  (s) 

Simull 

34 

0.0175 

139 

0.01 

20 

0.0075 

SimuT2 

18 

0.03 

58 

0.025 

14 

0.0075 

SimuB 

27 

0.025 

99 

0.01 

9 

0.0075 

Table  I.  Efficiency  on  a  single-processor  machine. 
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4.3.  RESULTS  ON  PARALLEL  MACHINES 

To  find  out  the  effect  of  parallelization  on  the  three  formulations  under  study,  the  stmulmton  Stmul  1 
has  been  executed  on  a  SUN  Sparc  Station  20  with  4  processors  HyperSparc  (RT625)  C  IOC 
with  256  Kb  of  cache  memoty.  Table  D  shows  the  reduction  in  CPU  lime  experienced  by  each 

method  as  the  number  of  processors  increases. 


#  of 

processors 

Augmented 

Lagrangian 

Classical 

Lagrangian 

Fully- 

Recursive 

1 

0.0 

0.0 

0.0 

2 

22.7 

6.5 

0.0 

3 

29.6 

9.0 

0.0 

4 

29.6 

9.0 

0.0 

Table  II.  Reduction  of  CPU  time  in  %  due  to  parallelization. 


It  may  be  seen  from  Table  D  that  the  method  that  takes  most  advantage  of  the  parallelization,  up  to 
30%  savings,  is  the  augmented  Lagrangian,  since  it  uses  it  during  the  matrix  formation  and  solution. 
The  classical  Lagrangian  only  takes  advantage  of  the  parallelization  at  the  time  of  solving  the 
equations,  achieving  up  to  9%  reduction  in  CPU  time.  On  the  other  hand,  the  fully-recursive 
formulation  does  not  profit  at  all  from  the  parallel  processing,  thus  indicating  to  be  a  method  suitable 
for  sequential  processing.  It  must  be  noticed  that,  in  all  cases,  the  fourth  processor  does  not 

contribute  to  improve  the  CPU  time. 


4.4  INFLUENCE  OF  NUMERICAL  STIFFNESS 

The  simulation  Simull  has  been  carried  out  again  using  the  three  formulations.  However,  this  time  it 
has  been  called  SimullS,  since  the  rubber  bumpers  of  stiffness  107  N/m  have  been  mtroduced  into 
the  four  suspensions  of  the  model,  so  that  the  algorithm  behavior  under  conditions  of  high  numerical 
stiffness  may  be  observed.  Table  III  compares  the  results  obtained  by  the  three  methods  when 
performing  Simull  and  SimullS. 


Augmented 

Lagrangian 

Classical  Lagrangian 

Fully-R 

ecursive 

CPU  (s) 

At  (s) 

CPU  (s) 

At  (s) 

CPU  (s) 

At  (s) 

Simull 

34 

0.0175 

139 

0.01 

20 

0.0075 

SimullS  . 

78 

0.005 

43  S 

0.0025 

3  4 

0.0025 

Table  III.  Influence, of  numerical  stiffness. 


Table  III  shows  that  numerical  stiffness  makes  all  the  methods  reduce  the  time  step  size.  The 
fully-recursive  by  a  factor  of  3,  the  augmented  Lagrangian  by  a  factor  of  3.5,  and  the  classical 
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Lagrangian  by  a  factor  of  4.  Thus  Ore  last  two  lose  penance  with  respect  to  the  fuUy^u  stve 

formulation.  The  penally  factor  of  the  augmented  Lagrangian  procedure  needs  to  be  tncrea 

,0  achieve  convergence  and,  at  lower  time  steps,  as  for  tnstance  10.  seconds,  thts  forrnu  a  on  s  r  s 

suffering  from  ill-conditioning  effects.  When  numerical  stiffness  remains  low,  thrs  problem  arises  ata 
1  &  _  -  .  ■  ■  .  • _ i  -iiffnarr  m oVpc  it  annppir  at 


step  size  of  10"4  or  1(T5  seconds,  but  the  presence  of  high  numerical  stiffness  makes  it  appear  at 


larger  time  steps. 


5.  Conclusions 

In  thts  paper  a  thorough  and  novel  companson  of  different  leadmg  methods  of  dynamic  simulation 

has  been  earned  out  for  large  multibody  systems.  The  conclustons  that  have  been  reached  may  be 

summarized  as  follows: 

.  -the  modeling  process  required  by  the  fuUy-recursive  method  leads  to  a  relative  coordmate 
formulation  that  ts  suitable  for  open  chain  configurations,  but  presents  serious  complexmes  and 
difficulties  tn  the  case  of  closed  kinematic  chains.  The  modeling  for  the  Lagrang.an  formulations  is 
much  simpler  and  may  be  done  with  dependent  coordinates.  Therefore  the  former  calls  for  special 
purpose  implementations  whereas  the  latter  is  more  suitable  for  general  purpose  simulators.  In 
respect  to  modeling,  an  analogy  may  be  established  between  the  dependent  coordinates  and  the 
finite  element  formulation  for  structural  analysis. 

•  The  augmented  Lagrangian  method  presents  a  leading  matrix  that  is  symmetric  and  positive 
definite,  thus  making  it  very  robust  to  address  simulations  with  singular  configurations.  Also,  this 
method  permits  the  differentiation  of  the  dynamic  equations  to  get  the  tangent  matrix  in  a  Newton- 
Raphson  iteration,  which  allows  for  convergence  at  large  time  step  sizes.  Furthermore,  it  takes  a 
serious  advantage  of  parallel  computing.  However,  it  suffers  from  ill-conditioning  when  the  time 
step  size  becomes  small,  thus  having  convergence  difficulties  for  stiff  systems.  A  precision 
environment  of  32-digits  could  be  the  cure  for  this  numerical  drawback. 

•  The  classical  Lagrangian  formulation  generates  an  almost  double-size  problem  than  its  augmented 
counterpart,  and  takes  less  profit  of  parallel  computing,  thus  performing  notably  worse.  In 
addition,  its  leadmg  matrix  is  less  robust,  causing  the  failure  of  the  integration  process  at  singular 
configurations  and  with  redundant  constraints. 

•  The  fuUy-recursive  procedure  performs  weU  for  stiff  and  non-stiff  systems.  Its  involved 
formulation  only  permit  a  fixed-point  iteration  scheme,  so  that  the  time  step  size  must  be  kept  small 
to  achieve  convergence.  Despite  of  this  fact,  it  shows  to  be  the  most  efficient  method  among  the 
three  compared,  but  does  not  obtain  any  advantage  from  parallel  computing. 

•  According  to  the  results  presented- in  this  paper,  real-time  simulation  of  large  multibody  systems 
can  be  achieved  with  medium-size  workstations,  opening  a  wide  field  for  development  of 
hardware-in-the-loop  and  man-in-the-loop  facilities,  virtual  prototyping  tools  and  intelligent  vehicle 

control  systems. 
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